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5.  Introduction 

Widely  recognized  strengths  of  ultrasonic  imaging  techniques  for  diagnosis  and 
monitoring  for  recurrence  of  breast  cancer  are  the  nonionizing  nature  of  acoustic 
waves,  the  ability  to  provide  good  contrast  between  fluid  and  parenchymal  tissue, 
portability,  and  good  patient  tolerance.  Ultrasonic  imaging  is  inexpensive  compared  to 
magnetic  resonance  imaging  and  is  competitive  in  cost  with  x-ray  mammography. 
However,  these  strengths  are  presently  offset  by  significant  weaknesses  that  include 
low  resolution,  inability  to  distinguish  between  solid  masses,  and  the  lack  of  quantitative 
information  about  breast  tissue  characteristics. 

For  some  years,  recognition  has  existed  that  high-resolution  and  quantitative  ultrasonic 
imaging  has  the  potential  to  provide  diagnostic  information  not  available  through 
conventional  ultrasonic  techniques  such  as  b-scanning.  Aberration  correction  promises 
to  increase  significantly  the  resolution  of  ultrasonic  mammography  so  that  small  lesions, 
tumors,  and  microcalcifications  can  be  clearly  imaged  even  through  fatty  or  dense 
breast  tissue.  Quantitative,  high-resolution  ultrasonic  imaging  is  particularly  desirable 
for  the  diagnosis  of  breast  disease  since  the  approach  has  the  potential  to  differentiate 
between  types  of  cancerous  tissue,  to  detect  and  characterize  microcalcifications,  to 
detect  small  tumors  even  in  dense  breast  tissue,  and  to  yield  diagnostically  useful 
images  without  operator-dependent  adjustments  of  parameters. 

The  purpose  of  this  research  is  to  develop  new  ultrasonic  imaging  concepts  based  on 
significant  advances  in  adaptive  focusing  and  breakthroughs  in  inverse  scattering  as 
well  as  on  newly  available  novel  experimental  apparatus.  To  focus  adaptively,  transmit 
signals  are  synthesized  from  scattered  signals  by  using  backpropagation  followed  by 
time-shift  estimation  and  compensation.  To  reconstruct  quantitatively,  eigenfunctions 
and  eigenvalues  of  an  operator  associated  with  scattering  measurements  are  used  to 
obtain  efficiently  images  of  distributed  scattering  objects. 
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6.  Body 

a.  Method 

The  method  employed  to  achieve  the  goal  of  quantitative  and  high-resolution  ultrasonic 
imaging  employs  a  combination  of  hardware  innovations,  careful  measurements  of 
ultrasonic  aberrations  and  scattering,  and  implementation  of  new  ideas  for  aberration 
correction  as  well  as  for  image  reconstruction  via  inverse  scattering.  The  novel 
hardware  consists  of  two  different  systems.  One  is  a  custom  made  pulse-echo 
apparatus  in  which  the  size  of  the  transmit  aperture  can  easily  be  altered  and  in  which 
scattering  objects  can  be  easily  changed  without  disturbing  the  aberrator.  The  other 
apparatus  consists  of  a  ring  transducer  array  with  2048  elements  and  associated 
circuitry  for  multiplexing  and  signal  control.  Measurements  performed  with  this  ring 
produce  experimental  data  from  which  the  degrading  effects  of  wavefront  distortion 
imaging  were  studied.  Measurements  also  provided  experimental  data  from  phantoms 
with  known  acoustic  parameters  for  image  reconstruction.  Receive  beamforming  used 
backpropagation  and  other  methods  of  space-time  processing.  Measured  scattering 
data  with  three  hundred  and  sixty  degrees  of  spatial  coverage  was  used  for  the 
reconstruction  of  quantitative  ultrasonic  images  through  the  use  of  newly  developed 
inverse  scattering  methods. 

b.  Results 

Ultrasonic  wavefront  distortion  produced  by  transmission  through  breast  tissue 
specimens  was  measured  in  a  two-dimensional  aperture.  Differences  in  arrival  time 
and  energy  level  between  the  measured  waveforms  and  references  that  account  for 
geometric  delay  and  spreading  were  calculated  as  was  the  waveform  similarity  factor. 
Details  of  this  investigation  were  reported  in  the  Journal  of  the  Acoustical  Society  of 
America.  A  copy  of  this  paper  is  included  as  Appendix  A.  In  this  paper,  maps  of 
distortion  are  shown  in  Figure  2  focussing  results  are  shown  in  Figure  3  and 
comparisons  with  corresponding  abdominal  wall  data  are  given  in  Figures  7  and  8. 

A  new  method  of  waveform  design  was  developed  and  implemented  to  produce  a 
spatially  limited  plane  wave  for  illumination  of  the  scattering  object  to  be  imaged.  The 
method  employed  backpropagation  of  the  desired  spatially  limited  plane  wavefront  to 
obtain  excitation  signals  for  the  ring  transducer  elements.  A  paper  describing  the 
methods  in  detail  was  published  in  the  IEEE  Transactions  on  Ultrasonics, 
Ferroelectrics,  and  Frequency  Control.  A  copy  of  this  paper  is  included  as  Appendix  B. 
Figure  6  in  this  paper  shows  illustrative  results. 

Estimation  and  correction  of  ultrasonic  wavefront  distortion  using  pulse-echo  data  has 
also  been  implemented  in  an  experimental  setting.  The  processing  used  random 
scattering  measured  in  a  specially  designed  pulse-echo  configuration  in  which  the 
transmit  beam  is  unperturbed  by  the  aberrating  medium.  Major  results  in  are: 
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1)  A  tightly  focussed  transmit  beam,  such  as  that  achieved  with  f/1.5,  is  needed  to 
produce  a  sufficiently  coherent  scattered  wavefront  so  that  time  delay  can  be 
accurately  estimated. 

2)  Accommodation  of  waveform  shape  distortion  in  time-delay  estimation,  such  as  in 
the  modified  least-mean-square-error  method  that  incorporates  a  global  smoothness 
requirement,  is  necessary  for  estimation  of  the  time-delay  surface  of  a  wavefront 
that  has  propagated  through  distributed  inhomogeneities. 

3)  The  isoplanatic  path  size  usually  decreases  with  increasing  severity  of  distortion  and 
is  influenced  by  compensation  technique. 

4)  Incorporation  of  waveform  shape  compensation,  such  as  by  using  backpropagation 
prior  to  time-shift  compensation,  improves  the  contrast  ratio  over  that  obtained  by 
time-shift  compensation  alone. 

A  paper  that  details  this  work  has  been  accepted  for  publication  in  the  IEEE 
Transactions  on  Ultrasonics,  Ferroelectrics,  and  Frequency  Control  and  a  preprint  of 
this  paper  is  included  as  Appendix  C.  Figures  9-14  of  the  preprint  show  results  that  are 
representative. 

A  novel  inverse  scattering  method  that  uses  eigenfunctions  of  the  scattering  operator 
was  examined  theoretically  and  using  computations.  A  unified  framework  that 
encompasses  eigenfunction  methods  of  focussing  and  image  reconstruction  in  arbitrary 
media  was  developed.  A  paper  detailing  this  work  has  been  published  in  the  Journal  of 
the  Acoustical  Society  of  America  and  is  included  as  Appendix  D  to  provide  more  detail. 
Figures  2  through  8  in  this  paper  show  representative  images. 

Representative  high-resolution  b-scan  images  produced  by  our  ring  transducer  system 
are  presented  in  Figure  1  along  with  a  corresponding  b-scan  image  produced  by  a 
conventional  ultrasonic  scanner  and  a  corresponding  x-ray  image  produced  by 
computed  tomography.  The  results  show  that  our  ring  transducer  system  is  capable  of 
producing  images  that  are  superior  to  those  of  commercial  b-scan  ultrasonic  imagers 
and  comparable  to  those  obtained  from  computed  x-ray  tomography. 
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Fig.  1.  Images  of  an  Anthropomorphic 
Breast  Phantom,  The  phantom  was 
imaged  in  a  plane  containing  two  9  mm 
diameter  cysts  on  the  periphery  of  the 
phantom  and  one  4  mm  diameter  cyst 
near  the  center.  Top  Left:  single  b-scan 
image  obtained  with  the  ring  transducer 
system  using  f/1.0  aperture  and  scan 
plane  diameter  of  80  mm.  Top  right:  b- 
scan  image  produced  by  a  state-of-the-art 
clinical  b-scanner  operating  at  the  same 
center  frequency  as  the  ring  transducer 
system.  Bottom  Left:  compounded  b- 
scan  image  from  six  geometrically 
corrected  views  with  the  ring  transducer 
system.  Bottom  Right:  X-ray  computed 
tomograph  obtained  with  a  state-of-the- 
art  scanner  using  a  1  mm  slice  thickness. 
These  images  show  that  ring  transducer 
b-scan  resolution  is  superior  to  current 
clinical  scanners  at  the  same  frequency 
and  is  comparable  with  x-ray  computed 
tomography. 


c.  Discussion 

The  studies  summarized  above  and  detailed  in  appendices  show  that  all  of  the  tasks  in 
the  statement  of  work  have  been  largely  accomplished  although  studies  employing 
breast  tissue  and  comparative  evaluations  of  imaging  modalities  were  more  limited  than 
originally  planned.  These  limitations  resulted  from  lost  time  caused  by  unforeseen 
equipment  failures  and  the  need  for  more  software  development  than  made  possible  by 
the  budget.  Nevertheless,  the  results  demonstrate  the  new  ultrasonic  imaging  methods 
developed  in  this  project  have  potential. 
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7.  Conclusions 


The  forgoing  results  and  discussion  show  that  the  research  generally  progressed  along 
the  lines  described  in  the  original  proposal.  Unanticipated  problems  limited  the  depth  of 
some  studies  and  the  extent  of  progress  but  important  theoretical,  experimental,  and 
computational  results  have  been  obtained.  Additional  studies  that  involve  more  use  of 
tissue  and  further  comparative  evaluation  of  imaging  modalities  are  needed  to 
determine  the  ultimate  capability  of  the  imaging  techniques  investigated  in  this  project 
to  improve  detection,  diagnosis,  and  treatment  monitoring  of  breast  cancer. 
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Ultrasonic  wavefront  distortion  produced  by  transmission  through  breast  tissue  specimens  was 
measured  in  a  two-dimensional  aperture.  Differences  in  arrival  time  and  energy  level  between  the 
measured  waveforms  and  references  that  account  for  geometric  delay  and  spreading  were 
calculated.  Also  calculated  was  a  waveform  similarity  factor  that  is  decreased  from  1.0  by  changes 
in  waveform  shape.  For  nine  different  breast  specimens,  the  arrival  time  fluctuations  had  an  average 
(±s.d.)  rms  value  of  66.8  (±12.6)  ns  and  an  associated  correlation  length  of  4.3  (±1.1)  mm,  while 
the  energy  level  fluctuations  had  an  average  rms  value  of  5.0  (±0.5)  dB  and  a  correlation  length  of 
3.4  (±0.8)  mm.  The  corresponding  waveform  similarity  factor  was  0.910  (±0.023).  The  effect  of 
the  wavefront  distortion  on  focusing  and  the  ability  of  time-shift  compensation  to  remove  the 
distortion  were  evaluated  by  comparing  parameters  such  as  the  —  30-dB  effective  radius,  the 
—  10-dB  peripheral  energy  ratio,  and  the  level  at  which  the  effective  radius  departs  from  an  ideal  by 
10%  for  the  focus  obtained  without  compensation,  with  time-shift  estimation  and  compensation  in 
the  aperture,  and  with  time-shift  estimation  and  compensation  performed  after  backpropagation.  For 
the  nine  specimens,  the  average  —10-dB  peripheral  energy  ratio  of  the  focused  beams  fell  from  3.82 
(±1.83)  for  the  uncompensated  data  to  0.96  (±0.18)  with  time-shift  compensation  in  the  aperture 
and  to  0.63  (±0.07)  with  time-shift  compensation  after  backpropagation.  The  average  -30-dB 
effective  radius  and  average  10%  deviation  level  were  4.5  (±0.8)  mm  and  —19.2  (±3.5)  dB, 
respectively,  for  compensation  in  the  aperture  and  3.2  (±0.7)  mm  and  —22.8  (±2.8)  dB, 
respectively,  for  compensation  after  backpropagation.  The  corresponding  radius  for  the 
uncompensated  data  was  not  meaningful  because  the  dynamic  range  of  the  focus  was  generally  less 
than  30  dB  in  the  elevation  direction,  while  the  average  10%  deviation  level  for  the  uncompensated 
data  was  —4.9  (±4.1)  dB.  The  results  indicate  that  wavefront  distortion  produced  by  breast 
significantly  degrades  ultrasonic  focus  in  the  low  MHz  frequency  range  and  that  much  of  this 
degradation  can  be  eliminated  using  wavefront  backpropagation  and  time-shift  compensation. 

PACS  numbers:  43.80.Cs,  43.80.Ev,  43.80.Vj 


INTRODUCTION 

Widely  recognized  strengths  of  ultrasonic  imaging  tech¬ 
niques  for  diagnosis  and  monitoring  of  breast  disease  are  the 
nonionizing  nature  of  acoustic  waves  and  the  ability  to  pro¬ 
vide  good  contrast  between  fluids  and  parenchymal  tissues. 
However,  despite  advances  in  transducer  technology,  breast 
ultrasonography  has  thus  far  been  relegated  to  ancillary  use, 
largely  because  resolution  is  inadequate.^’^  This  has  led  to 
consideration  of  the  limitations  imposed  on  ultrasonic  imag¬ 
ing  of  the  breast  by  wavefront  distortion  that  arises  from 
propagation  through  breast  tissue. 

Several  researchers  have  studied  ultrasonic  distortion 
produced  by  the  breast  and  considered  how  this  distortion 
may  be  compensated.  An  early  in  vivo  pulse-echo  study  by 
Moshfeghi  and  Waag^  showed  that  increasing  the  aperture  to 
//l.O  from  //2.6  for  breast  produced  only  about  f  the  resolu¬ 
tion  improvement  predicted  in  a  homogeneous  medium.  Tra- 
hey  et  al^  made  one-dimensional  transmission  measure¬ 


ments  of  ultrasonic  phase  distortion  caused  in  vivo  by 
propagation  through  breast  and  found  an  average  rms  arrival 
time  aberration  of  36.1  ns  for  22  subjects.  The  same  group 
later  extended  their  transmission  measurements  to  two  di¬ 
mensions  and  obtained  an  rms  value  of  55.3  ns  for  seven 
volunteers.^  They  concluded  that  phase  distortion  should  be 
measured  and  corrected  in  two  dimensions  but  did  not  men¬ 
tion  amplitude  distortion.  In  other  in  vivo  one-dimensional 
transmission  measurements,  Zhu  and  Steinberg  observed  se¬ 
vere  amplitude  distortion,  which  they  attributed  to  refraction 
in  addition  to  scattering, and  obtained  a  relation  between 
the  average  sidelobe  floor  and  the  normalized  variance  of  the 
amplitude  distortion  produced  by  the  breast.®’^  They  con¬ 
cluded  that  large  two-dimensional  arrays  and  new  algorithms 
that  correct  both  phase  and  amplitude  distortion  in  two  di¬ 
mensions  may  be  needed  to  reduce  the  effects  of  distortion 
produced  by  the  breast. 

Two  basic  algorithms  for  correction  of  ultrasonic  distor- 
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tion  have  been  investigated  for  general  imaging  applications. 
The  first  relies  on  cross  correlation  of  signals  in  the  aperture 
for  estimation  of  pulse  arrival  time/^’^^  while  the  other  ad¬ 
justs  beamformer  delays  to  maximize  signal  brightness. 
However,  neither  of  these  methods  addresses  the  problems  of 
amplitude  and  waveform  distortion,  since  both  use  only 
time-shift  compensation  in  the  aperture.  Another  technique 
has  been  investigated  for  the  removal  of  amplitude  and 
waveform  distortion  as  well  as  time-shift  distortion  in  spe¬ 
cialized  applications  when  a  point  source  is  present  as  may 
be  the  case  in  lithotripsy.  This  technique  employs  time- 
reversed  signals  but  is  limited  because  a  suitable  point  source 
is  not  generally  available  in  every  isoplanatic  patch  to  be 
imaged.  A  more  recent  method  developed  by  Liu  and 
Waag^^  can  be  used  to  remove  time-shift,  amplitude,  and 
waveform  distortion  in  general  imaging  applications.  Their 
approach  models  the  distorting  medium  as  a  phase  screen 
placed  some  distance  from  the  receiving  aperture  and  re¬ 
moves  amplitude  and  waveform  distortion  by  backpropaga- 
tion  of  the  wavefront  before  applying  time-shift  compensa¬ 
tion. 

This  paper  reports  a  study  of  ultrasonic  wavefront  dis¬ 
tortion  produced  by  breast  and  the  effectiveness  of  the  back- 
propagation  method  in  removing  the  distortion.  In  the  study, 
wavefronts  perturbed  by  transmission  through  breast  tissue 
were  measured  in  a  two-dimensional  aperture.  Statistics  de¬ 
scribing  arrival  time  variations,  energy  level  fluctuations,  and 
wave  shape  distortion  were  calculated  and  compared  to  val¬ 
ues  from  analogous  measurements  using  abdominal  wall.^^ 
The  received  waveforms  were  then  focused  without  compen¬ 
sation,  with  time-shift  estimation  and  compensation  in  the 
measurement  aperture,  and  with  backpropagation  followed 
by  time-shift  estimation  and  compensation  to  determine  the 
effectiveness  of  time-shift  compensation  with  and  without 
wavefront  backpropagation  for  the  improvement  of  focusing. 

I.  METHOD 

Breast  tissue  specimens  were  obtained  fresh  from  reduc¬ 
tion  mammoplasty  surgery  and  were  stored  frozen  if  not  used 
immediately  for  measurements.  The  specimens  came  from 
regions  of  the  breast  away  from  the  nipple  and  consisted  of 
fat,  glandular  and  connective  tissue,  and  a  surface  covering 
of  skin.  Each  specimen  was  essentially  planar  and  had  a 
surface  area  of  at  least  7X11  cm^.  The  average  thickness  was 
26,9  mm.  The  tissue  donors  were  women  ranging  in  age 
from  18  to  65  with  a  mean  age  of  34  years. 

Measurements  were  carried  out  using  the  procedure  de¬ 
tailed  in  Ref.  17  and  summarized  here  for  convenience.  A 
breast  tissue  specimen  was  placed  in  the  specimen  holder 
with  the  skin  facing  the  direction  of  the  receiving  transducer 
and  pressurized  to  500  psi  for  30  min  in  order  to  dissolve  any 
gas  bubbles  present  in  the  tissue.  The  specimen  holder  was 
then  mounted  in  the  experimental  chamber,  which  was  main¬ 
tained  at  37±1  throughout  the  measurement.  Ultrasonic 
pulses  emitted  by  a  hemispheric  transducer  were  received  by 
a  128-element  linear  array  immediately  after  propagation 
through  the  specimen.  A  two-dimensional  area  was  scanned 
by  translating  the  array  in  the  elevation  direction  using  an 
automated  stage.  At  each  elevation,  the  array  elements  were 
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accessed  sequentially  by  a  multiplexer.  The  signal  from  each 
element  was  digitized  into  12-bit  samples  for  a  period  of  11.8 
yLts  at  a  rate  of  20  MHz.  The  signal  was  recorded  19  times  at 
each  element  to  permit  noise  reduction  through  signal  aver¬ 
aging.  The  nominal  center  frequency  was  3.75  MHz  for  each 
transducer  and  the  ~6-dB  bandwidth  of  the  received  pulse 
was  about  2.2  MHz.  The  element  pitch  in  the  receiving  trans¬ 
ducer  was  0.72  mm  and  a  reflecting  mask  reduced  the  receiv¬ 
ing  elevation  to  1.44  mm.  A  period  of  about  35  min  was 
required  to  record  the  signals  from  all  128  array  elements  at 
each  of  32  elevations  spaced  1.44  mm  apart  for  a  total  of 
4096  positions  over  a  92.16X46.08  mm^  aperture.  The  total 
source-receiver  separation  was  about  165  mm  and  the 
specimen-receiver  gap  was  about  8  mm. 

Three  groups  of  measurements  were  made.  The  primary 
set  consisted  of  independent  measurements  using  nine  differ¬ 
ent  specimens.  In  addition,  two  specimens  were  employed 
for  sequential  measurements  in  which  a  1-cm  layer  was  re¬ 
moved  from  the  bottom  of  the  specimen  before  each  subse¬ 
quent  measurement.  For  two  other  specimens,  a  pair  of  mea¬ 
surements  was  made  with  the  source  in  each  of  two  positions 
located  12  mm  apart  in  the  array  direction. 

Variations  in  pulse  arrival  time  were  calculated  using  the 
reference  waveform  method  for  arrival  time  estimation  and 
differences  in  shape  among  the  received  waveforms  were 
quantified  using  the  waveform  similarity  factor  as  described 
in  Ref.  16.  Energy  level  fluctuations  were  calculated  as  de¬ 
scribed  in  Ref.  17.  The  calculations  are  outlined  here  to  iden¬ 
tify  the  main  steps.  A  reference  waveform  was  constructed 
for  a  set  of  data  as  an  average  of  all  the  waveforms  meeting 
a  cross-correlation  criterion  for  similarity.  The  reference 
pulse  was  then  cross  correlated  with  each  of  the  original 
waveforms  and  an  arrival  time  surface  was  calculated  from 
the  peaks  of  the  correlation  functions.  Smoothing  was  done 
to  remove  questionable  outlying  points.  A  two-dimensional, 
fourth-order  polynomial  fit  to  the  estimated  arrival  times  was 
used  to  position  a  window  on  the  original  waveforms,  and 
arrival  time  estimation  was  repeated  using  the  windowed 
data.  Geometric  effects  were  removed  by  fitting  a  two- 
dimensional,  fourth-order  polynomial  to  the  newly  calculated 
arrival  time  surface  and  subtracting  the  result  to  obtain  the 
arrival  time  fluctuations.  Fluctuations  in  energy  level  across 
the  received  wavefront  were  calculated  by  summing  the 
squared  amplitudes  of  the  windowed  signals  at  each  position, 
converting  the  results  to  decibel  units,  and  subtracting  a  fit¬ 
ted  two-dimensional,  fourth-order  polynomial  from  the  re¬ 
sult,  This  computation,  however,  did  not  employ  the  20-dB 
restriction  on  the  dynamic  range  as  in  Ref,  17,  since  no  ex¬ 
ceptionally  low-energy  values  occurred.  The  waveform  simi¬ 
larity  factor,  which  ranges  from  an  ideal  value  of  1  to  a 
minimum  of  0  and  is  insensitive  to  absolute  amplitude  and 
arrival  time  like  a  correlation  coefficient,  was  computed  for 
all  the  windowed  waveforms  throughout  the  aperture. 

Two  forms  of  compensation  for  wavefront  distortion 
were  applied  to  compare  their  effect  on  focusing.  In  one, 
time-shift  estimation  and  compensation  were  performed  di¬ 
rectly  in  the  receiving  aperture.  In  the  other,  time-shift  esti¬ 
mation  and  compensation  were  performed  after  the  wave- 
front  had  been  backpropagated  to  the  point  of  maximum 
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FIG.  1.  Representative  waveforms  compensated  for  geometric  travel  time,  (a)  H2O.  (b)  BRS  10.  (c)  BRS  8a.  (d)  BRS  8b.  Each  column  of  panels  shows 
waveforms  at  sequential  increments  of  2.88  mm  in  elevation  across  the  first  half  of  the  46.08-mm  aperture.  At  each  elevation,  the  horizontal  coordinate  is  the 
array  direction  and  spans  a  distance  of  92.16  mm  in  0.72-mm  increments  while  the  vertical  coordinate  is  time  and  spans  an  interval  of  2.0  fis  in  0.05-//S 
increments.  Signal  amplitude  is  shown  linearly  on  a  gray  scale  with  the  maximum  signal  in  the  two-dimensional  aperture  for  each  measurement  represented 
by  white  and  the  corresponding  negative  value  by  black. 


waveform  similarity.  The  backpropagation  used  the  angular 
spectrum  method  as  described  in  Ref.  16.  The  reference 
waveform  method  was  employed  in  both  forms  for  time-shift 
estimation. 

The  efficacy  of  the  two  correction  procedures,  as  well  as 
the  effect  of  the  original  distortion,  were  evaluated  by  focus¬ 
ing  each  set  of  data  at  180  mm  via  a  Fourier  transform 
method  as  described  in  Ref.  18.  Each  focus  was  described  by 
the  “10,  “20,  and  —30  dB  effective  radii,  which  are  half  the 
cube  root  of  the  product  of  the  corresponding  effective 
widths  in  the  array,  elevation,  and  time  directions.  The  effec¬ 
tive  radius  several  dB  down  from  the  peak  is  a  useful  mea¬ 
sure  of  point  resolution  while  the  effective  radius  many  dB 
down  from  the  peak  is  a  measure  of  contrast  resolution.  The 
—  10-dB  peripheral  energy  ratio,  which  is  the  ratio  of  energy 
outside  an  ellipsoid  bounded  by  the  —10  dB  effective  widths 
to  the  energy  inside  that  ellipsoid,  was  also  used  to  describe 
each  focus  quantitatively.  Since  the  energy  outside  a  speci¬ 
fied  region  around  the  main  peak  of  the  focus  can  be  viewed 
as  an  integration  of  sidelobe  intensity,  the  peripheral  energy 
ratio  is  another  measure  of  contrast  resolution.  Additionally, 
each  focus  was  described  by  the  10%  deviation  level,  which 
is  the  level  at  which  the  effective  radius  becomes  10%  larger 
than  that  produced  by  ideal  waveforms  obtained  for  the  data 


set  by  replicating  its  average  time-shift  compensated  wave¬ 
form  throughout  the  aperture.  This  level  provides  a  measure 
of  the  degree  to  which  the  central  region  of  actual  focus  is 
ideal. 

II.  RESULTS 

Representative  sets  of  waveforms  recorded  after  trans¬ 
mission  through  breast  tissue  specimens  and  corrected  for 
geometric  arrival  time  are  shown  in  Fig.  1  with  correspond¬ 
ing  representative  water  path  waveforms  to  illustrate  the 
range  and  combination  of  distortion  levels  encountered  in 
this  study.  The  water  path  waveforms  in  (a)  show  minimal 
arrival  time  and  energy  level  fluctuation  and  nearly  ideal 
waveform  similarity.  The  tissue  path  waveforms  in  (b)  ex¬ 
hibit  low  arrival  time  and  energy  level  fluctuation  and  mod¬ 
erate  waveform  distortion.  The  waveforms  in  (c)  have  mod¬ 
erate  arrival  time  variation,  high  energy  level  fluctuation,  and 
low  wave  shape  distortion.  The  waveforms  in  (d)  show  high 
arrival  time  and  waveform  distortion  but  moderate  energy 
level  fluctuation. 

Arrival  time  fluctuations  and  energy  level  fluctuations 
produced  by  nine  different  breast  tissue  specimens  are  shown 
in  Fig.  2  and  statistics  of  these  data  as  well  as  the  waveform 
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FIG.  2.  Breast  tissue  path  arrival  time  and  energy  level  fluctuations  for  nine 
different  specimens,  (a)  BRS  7.  (b)  BRS  8a.  (c)  BRS  8b.  (d)  BRS  9a.  (e) 
BRS  9b.  (f)  BRS  10.  (g)  BRS  11.  (h)  BRS  13.  (i)  BRS  17.  In  the  left  panel 
of  each  pair,  arrival  time  difference  is  shown  on  a  linear  scale  with  a  maxi¬ 
mum  arrival  time  fluctuation  of  +150  ns  represented  by  white  and  a  mini¬ 
mum  arrival  time  fluctuation  of  - 150  ns  represented  by  black.  In  the  right 
panel  of  each  pair,  energy  level  fluctuations  are  shown  on  a  log  scale  with  a 
maximum  positive  excursion  of  + 10  dB  represented  by  white  and  a  maxi¬ 
mum  negative  excursion  of  — 10  dB  represented  by  black.  In  all  panels,  the 
horizontal  coordinate  is  the  array  direction  and  spans  a  distance  of  92.16 
mm  in  0.72-mm  increments  while  the  vertical  coordinate  corresponds  to 
position  of  the  array  in  elevation  and  spans  a  distance  of  46.08  mm  with 
points  interpolated  from  measurements  at  1.44-mm  intervals  to  produce  data 
at  0.72-mm  increments. 
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similarity  factor  are  given  in  Table  I  for  each  set  of  measured 
waveforms.  The  arrival  time  and  energy  level  fluctuation  pat¬ 
terns  produced  by  a  given  specimen  are  roughly  similar  and 
the  specimen- to-specimen  variability  is  not  great.  Consistent 
with  these  observations,  the  effective  arrival  time  and  energy 
level  fluctuation  correlation  lengths  are  similar  for  each 
specimen  while  the  average  waveform  similarity  factor  has  a 
small  standard  deviation. 

Focal  plane  amplitude  at  representative  instants  of  time 
are  shown  in  Fig.  3  for  typical  measured  waveforms  that 
have  been  focused  without  compensation,  with  time-shift 
compensation  in  the  aperture,  and  with  time-shift  compensa¬ 
tion  following  backpropagation  along  with  the  focus  ob¬ 
tained  with  ideal  data.  The  spread  of  the  focus  obtained  from 
uncompensated  waveforms  is  large  relative  to  the  spread  of 
the  focus  obtained  from  ideal  data.  The  focus  of  the  wave¬ 
forms  that  have  been  time-shift  compensated  is  more  con¬ 
centrated  than  that  obtained  from  uncompensated  waveforms 
and  further  concentration  is  obtained  in  the  focus  of  wave¬ 
forms  that  have  been  time-shift  compensated  following  back- 
propagation,  but  the  focus  of  each  is  not  as  concentrated  as 
the  focus  obtained  from  ideal  waveforms. 

Effective  radius  curves  of  the  focus  obtained  from  the 
same  typical  measured  waveforms  for  uncompensated,  time- 
shift  compensated,  backpropagated  and  time-shift  compen¬ 
sated,  and  ideal  cases  are  shown  in  Fig.  4.  The  effective 
radius  of  the  focus  obtained  without  compensation  departs 
from  the  ideal  case  by  10%  at  a  level  —5.2  dB  below  the 
peak.  Time-shift  compensation  in  the  aperture  reduces  the 
10%  deviation  level  to  -21.4  dB  while  time-shift  compen¬ 
sation  following  backpropagation  has  a  10%  deviation  level 
of  -25.8  dB  and  has  an  effective  radius  that  is  appreciably 
narrower  at  levels  30  to  40  dB  below  the  peak. 

The  focus  obtained  with  each  of  the  nine  different  sets 
of  data  as  well  as  the  improvement  in  focus  obtained  using 
time-shift  compensation  in  the  aperture  and  using  time-shift 
compensation  after  backpropagation  are  described  by  the  pa¬ 
rameters  in  Table  IL  The  data  show  that  time-shift  compen¬ 
sation  in  the  aperture  and  time-shift  compensation  following 
backpropagation  result  in  an  effective  radius  at  the  focus  that 
is  similar  at  the  -10-  and  -20-dB  levels  but  that  time-shift 
compensation  following  backpropagation  improves  the 
— 30-dB  effective  radius  substantially  over  that  obtained  with 
time-shift  compensation  in  the  receiving  aperture.  This  indi¬ 
cates  that  time-shift  compensation  with  or  without  back- 
propagation  yields  about  the  same  point  resolution  but  that 
time-shift  compensation  after  backpropagation  improves  the 
contrast  resolution  more  than  time-shift  compensation  in  the 
receiving  aperture  does.  The  —  10-dB  peripheral  energy  ratio 
obtained  with  time-shift  compensation  in  the  aperture  aver¬ 
ages  about  25%  of  that  without  compensation  while  the 
-10-dB  peripheral  energy  ratio  obtained  with  time-shift 
compensation  following  backpropagation  is  about  16%  of 
that  obtained  without  compensation.  These  peripheral  energy 
ratios  are  another  indication  that  the  sidelobe  level,  and 
therefore  the  contrast  ratio,  is  improved  more  by  backpropa¬ 
gation  followed  by  time-shift  compensation  than  by  time- 
shift  compensation  alone. 

The  arrival  time  and  energy  level  fluctuations  produced 
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TABLE  I.  Breast  tissue  path  statistics  of  wavefront  distortion  produced  by  nine  different  specimens. 


Specimen 

number 

Specimen 

thickness 

(mm) 

Arrival  time  fluctuations 
rms  99,5%  Effective 

value  value  corr.  len. 

(ns)  (ns)  (mm) 

Energy  level  fluctuations 
rms  99.5%  Effective 

value  value  corr.  len. 

(dB)  (dB)  (mm) 

Waveform 

similarity 

factor 

7 

15-25 

63.1 

196.2 

5.82 

4.86 

14.21 

3.20 

0.940 

8a 

30-35 

66.5 

222.0 

4.88 

6.08 

16.35 

4.85 

0.926 

8b 

30 

85.6 

252.0 

5.16 

5.02 

12.99 

3.41 

0.869 

9a 

40 

76.5 

225,7 

4.81 

5.65 

15.10 

4,25 

0.914 

9b 

35 

79.5 

276.6 

4.20 

4.98 

13.09 

3.04 

0.883 

10 

20-25 

59.5 

212,8 

3.38 

4.67 

12.75 

2.87 

0.908 

11 

20-25 

70.7 

252.2 

4.73 

4.77 

12.44 

2.75 

0.903 

13 

15-20 

49.3 

160.8 

3.36 

4.72 

13.44 

3.39 

0.930 

17 

20-25 

50.5 

162.3 

2.40 

4.53 

12.46 

2.46 

0.918 

mean 

26.9 

66.8 

217.8 

4.30 

5.03 

13.65 

3.36 

0.910 

s.d. 

7.7 

12.6 

39.9 

1.07 

0.51 

1.33 

0.76 

0.023 

by  different  thicknesses  of  tissue  for  two  specimens  are 
shown  in  Fig.  5  and  described  statistically  in  Table  III.  The 
arrival  time  fluctuations,  energy  level  fluctuations,  and  wave¬ 
form  distortion  increase  with  specimen  thickness  for  all  but 
one  measurement.  The  arrival  time  and  energy  level  fluctua¬ 
tion  patterns  for  the  different  thicknesses  of  the  same  speci¬ 
men  are  seen  to  vary  considerably.  The  focus  characteristics 
associated  with  these  patterns  are  given  in  Table  IV. 

The  arrival  time  and  energy  level  fluctuations  in  the  ap¬ 
erture  and  after  backpropagation  are  shown  in  Fig.  6  for  two 
specimens  for  each  of  two  source  positions  12  mm  apart.  The 
features  in  the  fluctuation  patterns  obtained  for  each  speci¬ 
men  are  recognizable  but  shifted  in  the  patterns  obtained 
after  the  change  in  source  position.  The  features  are  also 
changed  but  still  identifiable  after  backpropagation.  The 
similarity  of  the  patterns  is  indicated  quantitatively  by  the 
distortion  statistics  given  in  Table  V  and  by  the  focus  char¬ 
acteristics  presented  in  Table  VI  for  compensation  using  the 
time  shifts  calculated  from  the  waveforms  produced  with  the 
source  at  the  focal  position.  The  data  in  Table  VI  also  show 
that  compensation  using  the  time  shifts  calculated  from 
waveforms  produced  with  the  source  12  mm  from  the  focal 
position  is  less  effective  than  compensation  using  the  time 
shifts  from  waveforms  produced  with  the  source  at  the  posi¬ 
tion  of  focus. 

III.  DISCUSSION 

The  specimens  studied  here  were  sections  of  breast  tis¬ 
sue  rather  than  whole  breasts.  They  came  from  unusually 
large  breasts  and  ranged  from  15-40  mm  in  thickness  with 
nearly  planar  surfaces.  Thus  while  the  measurements  illus¬ 
trate  elements  of  distortion  produced  by  the  breast,  they  may 
not  describe  conditions  typically  encountered  in  the  clinic. 
Nevertheless,  the  information  obtained  in  this  study  under 
precisely  controlled  experimental  conditions  adds  to  the 
available  data  about  the  ultrasonic  wavefront  distortion  pro¬ 
duced  by  the  breast  and  should  be  useful  in  the  development 
of  imaging  instruments  with  better  resolution  for  improved 
ultrasonic  diagnosis  of  breast  disease. 

Refraction  at  the  specimen-water  interfaces  could  con¬ 
tribute  to  the  distortion  measured  in  this  transmission  study. 
Prior  investigations^^’^^  have  shown  that  this  contribution  is 
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negligible  for  abdominal  wall,  for  which  the  surfaces  are  the 
skin  and  the  peritoneum.  The  breast  tissue  specimens,  how¬ 
ever,  had  a  cut  surface  on  their  lower  side,  so  a  special  ex¬ 
periment  was  undertaken  to  investigate  the  contribution  of 
refraction  at  this  surface  to  the  measured  distortion.  A  thick 
specimen  (BRS  9b)  was  selected  and  two  measurements 
were  made  over  approximately  the  same  area,  one  with  the 
skin  side  up  and  the  other  with  the  skin  side  down.  Because 
of  the  30-mm  difference  in  propagation  distance  from  the  cut 
surface  to  the  receiver  in  these  measurements,  differences  in 
the  distortion  patterns  from  these  measurements  should  be 
evident  if  refraction  at  the  cut  surface  were  an  appreciable 
source  of  the  distortion.  In  this  experiment,  however,  the 
arrival  time  fluctuation  rms  values  changed  by  less  than  9% 
while  the  energy  level  fluctuation  rms  values  and  all  corre¬ 
lation  length  values  differed  by  less  than  3%.  These  changes 
are  within  the  range  expected  from  reproducibility  studies^^ 
given  the  uncertainty  in  the  scan  position  and  the  differences 
in  propagation  path  caused  by  inverting  the  specimen.  The 
arrival  time  and  energy  level  fluctuation  maps  from  the  two 
measurements  were  also  remarkably  similar.  Therefore,  the 
cut  specimen  surface,  which  is  smoothed  by  the  kapton 
membrane  during  measurements,  is  not  considered  to  con¬ 
tribute  significantly  to  the  distortion  measured  for  the  breast 
specimens. 

The  data  processing  used  to  determine  the  wavefront 
distortion  in  this  study  employed  the  same  calculated  fourth- 
order  polynomial  references  for  the  determination  of  arrival 
and  energy  fluctuation  as  were  described  in  the  report^^  of 
wavefront  distortion  by  abdominal  wall,  but  incorporated 
other  improvements  and  extensions.  The  polynomial  fits,  as 
noted  in  Ref.  17,  avoid  the  need  for  accurate  knowledge 
about  tissue  dimensions  and  local  acoustic  characteristics 
and  compensate  for  small  misalignments  and  low  spatial- 
frequency  variations.  The  adaptive  reference  waveform 
method  used  to  estimate  arrival  time  differences  in  this  paper 
is  considered  more  accurate  than  the  least-mean-square  error 
method  employed  previously  because,  on  average,  fewer 
than  20%  of  the  cross  correlations  computed  for  each  ab¬ 
dominal  wall  data  set  using  the  adaptive  reference  waveform 
method  had  a  peak  value  below  0.8,  while  about  30%  of  the 
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FIG.  3.  Focal  plane  time  histories  of  representative  measured  data  (BRS  7).  Each  panel  shows  the  bipolar  distribution  of  signal  amplitude  as  a  shade  of  gray 
on  a  50-dB  log  scale  for  each  polarity  in  the  focal  plane  at  an  instant  of  time.  In  all  the  panels,  the  horizontal  coordinate  is  elevation  and  spans  37.504  mm 
while  the  vertical  coordinate  is  in  the  array  direction  and  spans  56.256  mm.  The  number  beneath  each  panel  identifies  the  (zero-origin)  instant  of  time  in  the 
128-point  interval  employed  in  the  temporal  Fourier  transform.  First  row:  Uncompensated  data.  Second  row:  After  time-shift  estimation  and  compensation  in 
the  aperture.  Third  row:  After  backpropagation  of  40  mm  followed  by  time-shift  estimation  and  compensation.  Fourth  row:  Ideal  data. 


peaks  were  below  0.8  for  the  least-mean-square  error 
method.  Also,  the  reference  waveform  method  gives  consis¬ 
tently  lower  estimates  of  distortion  for  water  paths  and  yields 
superior  focus  characteristics  when  employed  for  time-shift 
compensation.  Therefore,  the  reference  waveform  method 
was  used  in  this  study  and  the  data  in  Ref.  16  were  repro¬ 
cessed  using  the  reference  waveform  method  to  permit  a 
comparison  of  the  breast  results  discussed  here  with  the  ab¬ 
dominal  wall  data. 

In  the  nine  independent  measurements  reported  here,  an 
average  (±s.d.)  of  2807  (±228)  waveforms  or  68.5%  of  the 
4096  received  waveforms  were  similar  enough  to  be  incor¬ 
porated  into  the  final  reference  waveform.  The  average  of  the 
maximum  value  for  correlations  between  each  of  the  win¬ 
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dowed  waveforms  and  the  final  reference  waveform  was 
0.868  (±0.018).  This  figure,  calculated  using  all  of  the  wave¬ 
forms  in  the  aperture  before  any  smoothing,  is  similar  to  the 
average  cross-correlation  value  of  0.86  calculated  for  neigh¬ 
boring  waveforms  in  the  35  data  sets  selected  for  analysis  by 
Freiburger  et  al^  in  their  two-dimensional  study  of  the 
breast.  An  average  of  758  (±183)  or  18.5%  of  the  calculated 
delays  were  deemed  unacceptable  because  they  deviated  sig¬ 
nificantly  from  the  overall  delay  surface.  Of  these,  495 
(about  two  thirds)  were  replaced  by  delays  from  new  peak 
searches  while  the  remaining  263  (about  one  third)  were  cor¬ 
rected  by  smoothing. 

The  representative  waveforms  in  Fig.  1  illustrate  the 
variability  of  the  wavefronts  from  which  the  characteristics 
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(FIG.  4.  Effective  radius  of  the  focus  obtained  with  representative  measured 
waveform  data  (BRS  7).  Unc= uncompensated  data.  TSC= time-shift  com¬ 
pensation  in  the  receiving  aperture.  BP  &  TSC=backpropagation  followed 
by  time-shift  compensation.  Ideal = replication  of  a  single  average  waveform 
throughout  the  receiving  aperture. 

of  wavefront  distortion  have  been  found.  The  uniformity  of 
the  water  path  waveforms  indicates  that  the  error  introduced 
by  the  apparatus  is  small  compared  to  the  fluctuations  pro¬ 
duced  by  the  tissue  path.  The  presence  of  tissue  creates 
waveform  shape  changes  which  decorrelate  the  waveforms 
and  introduce  uncertainty  in  the  estimation  of  arrival  time. 
Improvement  of  arrival  time  estimation  accuracy,  and  thus 
the  correction  of  wavefront  distortion,  necessitates  the  re¬ 
moval  of  wave  shape  distortion  and  has  motivated  the  intro¬ 
duction  of  a  backpropagation  step  before  time-shift  compen¬ 
sation. 

!  The  arrival  time  and  energy  level  fluctuation  maps  in 

Fig.  2  for  the  nine  independent  breast  tissue  measurements 
show  the  range  of  ultrasonic  wavefront  distortion  patterns 
,  produced  by  human  breast  tissue  in  this  study.  These  patterns 

have  different  characteristics  from  those  reported^^  for  the 
abdominal  wall.  For  example,  the  breast  arrival  time  fluctua¬ 
tion  maps  contain  smaller  features  than  those  in  the  abdomi- 
i  nal  wall  arrival  time  plots.  The  backgrounds  of  the  breast 

arrival  time  and  energy  level  fluctuation  maps  are  also  com¬ 
prised  of  larger,  more  irregular  patches  than  those  in  the 
abdominal  wall  plots.  However,  as  is  the  case  with  the  ab¬ 
dominal  wall  patterns,  some  breast  arrival  time  fluctuation 


FIG.  5.  Breast  tissue  path  arrival  time  and  energy  level  fluctuations  for 
reductions  of  specimen  thickness,  (a)-(c)  BRS  9a.  (d)  and  (e)  BRS  11.  For 
each  specimen,  the  results  are  presented  from  top  to  bottom  in  order  of 
decreasing  thickness.  The  format  and  scales  in  each  pair  of  panels  are  the 
same  as  in  Fig.  2. 

map  features  appear  to  correlate  with  energy  level  fluctua¬ 
tions  while  others  do  not. 

The  bar  charts  in  Fig.  7  compare  the  arrival  time,  energy 
level,  and  waveform  distortion  statistics  for  the  nine  breast 
specimen  measurements  to  those  in  Ref.  17  for  abdominal 
wall  after  recalculation  using  the  reference  waveform 
method.  In  general,  breast  tissue  appears  to  cause  more  dis¬ 
tortion  than  the  abdominal  wall.  The  rms  arrival  time  fluc¬ 
tuations  produced  by  the  breast  specimens  in  this  study  have 


TABLE  II.  Breast  tissue  path  focus  characteristics  for  nine  different  specimens,  effective  radius.  PER = peripheral  energy  ratio.  Unc=uncompensated. 
TSC=time-shift  compensation.  BP = backpropagation  followed  by  TSC.  BP  dist==  distance  of  backpropagation  for  maximum  waveform  similarity. 


Specimen 

number 

BP  dist. 
(mm) 

Unc 

(mm) 

-10  dB 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

-20  dB  r, 
TSC 
(mm) 

BP 

(mm) 

Unc 

(mm) 

-30  dB 
TSC 
(mm) 

BP 

(mm) 

Unc 

10  dB  PER 

TSC  BP 

Unc 

(dB) 

10%  dev.  lev. 

TSC  BP 

(dB)  (dB) 

7 

30 

1.38 

1.04 

1.03 

3.94 

1.52 

1.49 

3.77 

2.60 

3.350 

0.766 

0.524 

-5.2 

-21.4 

-25.8 

8a 

40 

1.48 

1.09 

1.06 

4.39 

1.93 

1.69 

4.66 

3.44 

2.733 

0.905 

0.629 

-1.0 

-12.8 

-17.0 

8b 

50 

1.67 

1.07 

1.07 

7.60 

1.92 

1.54 

5.41 

3.24 

8.031 

1.280 

0.595 

-0.9 

-17.7 

-22.0 

9a 

50 

2.16 

1.06 

1.04 

5.02 

2.06 

1.55 

4.48 

3.76 

3.966 

1.128 

0.654 

-4.9 

-15.2 

-20.9 

9b 

30 

1.52 

1.12 

1.10 

5.73 

1.63 

1.57 

6.15 

4.49 

4.943 

1.044 

0.742 

-1.0 

-20.2 

-23.7 

10 

40 

1.17 

1.04 

1.03 

3.39 

1.59 

1.50 

3.64 

2.54 

3.957 

0.968 

0.592 

-8.0 

-19.4 

-22.2 

11 

30 

2.78 

1.04 

1.02 

5.70 

1.50 

1.47 

4.34 

3.03 

2.618 

0.962 

0.700 

-3.4 

-22.2 

-24.4 

13 

30 

1.16 

1.05 

1.03 

2.71 

1.54 

1.48 

5.92 

3.74 

2.89 

1.766 

0.712 

0.554 

-6.6 

-20.1 

-23.5 

17 

30 

1.04 

1.01 

1.00 

2.83 

1.43 

1.40 

7.52 

4.25 

2.57 

3.038 

0.877 

0.631 

-13.3 

-24.0 

-26.1 

mean 

1.60 

1.06 

1.04 

4.59 

1.68 

1.52 

4.49 

3.17 

3.822 

0.960 

0.625 

-4.9 

-19.2 

-22.8 

s.d. 

0.56 

0.03 

0.03 

1.59 

0.23 

0.08 

0.83 

0.65 

1.827 

0.176 

0.068 

4.1 

3.5 

2.78 
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TABLE  in.  Breast  tissue  path  statistics  of  wavefront  distortion  for  reductions  of  specimen  thickness. 


Arrival  time  fluctuations  Energy  level  fluctuations 


Specimen 

number 

Specimen 

thickness 

(mm) 

rms 

value 

(ns) 

99.5% 

value 

(ns) 

Effective 

corn  len. 
(mm) 

rms 

value 

(dB) 

99.5% 

value 

(dB) 

Effective 
corn  len. 
(mm) 

Waveform 

similarity 

factor 

9a 

40 

76.5 

225.7 

4.81 

5.65 

15.10 

4.25 

0.914 

30 

58.2 

183.6 

4.86 

5.44 

14.79 

4.06 

0.947 

18 

54.6 

157.6 

5.12 

5.01 

15.31 

2.92 

0.965 

11 

20-25 

70.7 

252.2 

4.73 

4.77 

12.44 

2.75 

0.903 

15 

56.2 

168.1 

4.10 

4.92 

14.67 

2.80 

0.933 

a  mean  of  66.8  ns,  which  is  about  20%  higher  than  the  mean 
of  55.5  ns  for  the  14  abdominal  wall  samples.  The  rms  en¬ 
ergy  level  fluctuations  produced  by  breast  in  this  study  have 
a  mean  value  of  5.03  dB,  which  is  about  39%  higher  than  the 
mean  of  3.62  dB  for  the  abdominal  wall.  Wave  shape  distor¬ 
tion  is  also  greater  for  the  breast  studies,  which  have  a  mean 
waveform  similarity  factor  of  0.910  versus  0.938  for  the  ab¬ 
dominal  wall. 

The  difference  in  the  average  thickness  of  breast  and 
abdominal  wall  specimens  may  contribute  to  the  difference 
in  the  average  arrival  time  fluctuation  rms  values  for  the  two 
types  of  specimens  but  appears  to  be  a  secondary  source  of 
the  variation  in  energy  level  fluctuations.  The  mean  thickness 
of  the  breast  specimens  is  25%  larger  than  that  of  the  ab¬ 
dominal  walls  and  the  average  arrival  time  fluctuation  of  the 
breast  specimens  is  about  20%  larger  than  that  of  the  ab¬ 
dominal  walls.  In  addition,  the  arrival  time  fluctuation  rms 


values  measured  for  specimens  of  both  types  having  similar 
thicknesses  were  approximately  the  same.  Thus  the  greater 
thickness  of  the  breast  specimens  likely  contributed  impor¬ 
tantly  to  the  greater  arrival  time  fluctuations  observed  in  this 
study.  However,  a  comparison  of  similarly  thick  breast  and 
abdominal  wall  specimens  revealed  that  the  energy  level 
fluctuation  rms  values  of  the  breast  specimens  were  larger  in 
every  case.  This  suggests  that  breast  tissue  produces  more 
energy  fluctuation  per  unit  thickness  than  does  abdominal 
wall  tissue,  so  that  the  mean  rms  energy  level  fluctuation 
produced  by  the  breast  specimens  would  be  greater  than  that 
of  the  abdominal  walls  even  without  the  difference  in  their 
thicknesses. 

The  spatial  distributions  of  arrival  time  and  energy  level 
fluctuations  caused  by  breast  and  abdominal  wall  tissue  dif¬ 
fer  as  well.  For  the  breast,  the  effective  correlation  length  of 
arrival  time  fluctuation  has  a  mean  of  4.31  mm,  which  is 


(a) 


(b) 


FIG.  6.  Arrival  time  and  energy  level  fluctuations  before  and  after  backpropagation  for  two  source  positions  12  mm  apart,  (a)  BRS  8a.  (b)  BRS  17.  For  each 
specimen,  the  four-panel  set  on  the  left  shows  arrival  time  fluctuations  and  the  four-panel  set  on  the  right  shows  the  corresponding  energy  level  fluctuations. 
In  each  four-panel  set,  the  upper  row  shows  the  fluctuations  with  the  source  to  the  right  and  the  lower  row  shows  the  fluctuations  with  the  source  to  the  left, 
while  the  left  column  shows  the  fluctuations  in  the  measurement  aperture  and  the  right  column  shows  the  corresponding  fluctuations  after  backpropagation  to 
the  distance  of  maximum  waveform  similarity.  The  scales  in  the  panels  are  the  same  as  in  Fig.  2. 
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TABLE  IV.  Breast  tissue  path  focus  characteristics  for  reductions  of  specimen  thickness.  r^=effective  radius.  PER = peripheral  energy  ratio.  Unc 
= uncompensated.  TSC= time-shift  compensation.  BP=backpropagation  followed  by  TSC.  BP  Dist=distance  of  backpropagation  for  maximum  waveform 
similarity. 


Specimen 

number 

BP  dist. 
(mm) 

Unc 

(mm) 

-10  dB 

TSC 

(mm) 

BP 

(mm) 

Unc 

(mm) 

-20  dB  r, 
TSC 
(mm) 

BP 

(mm) 

Unc 

(mm) 

-30  dB 

TSC 

(mm) 

BP 

(mm) 

Unc 

10  dB  PER 

TSC  BP 

Unc 

(dB) 

10%  dev.  lev. 

TSC  BP 

(dB)  (dB) 

9a 

50 

2.16 

1.06 

1.04 

5.02 

1.06 

1.55 

4.48 

3.76 

3.966 

1.128 

0.654 

-4.9 

-15.2 

-20.9 

40 

1.14 

1.04 

1.03 

3.66 

1.73 

1.56 

7.20 

3.82 

2.89 

2.372 

0.830 

0.540 

-1.0 

-17.0 

-19.8 

20 

1.56 

1.03 

1.03 

4.10 

1.49 

1.46 

8.06 

4.21 

2.94 

2.954 

0.673 

0.508 

-1.0 

-22.1 

-25.0 

11 

30 

2.78 

1.04 

1.02 

5.70 

1.50 

1.47 

4.34 

3.03 

2.618 

0.962 

0.700 

-3.4 

-22.2 

-24.4 

20 

1.36 

1.01 

1.01 

3.88 

1.43 

1.41 

2.99 

2.52 

2.779 

0.831 

0.565 

-5.6 

-22.2 

-24.2 

26%  smaller  than  the  5.78  mm  mean  calculated  for  the  ab¬ 
dominal  wall.  The  effective  correlation  length  of  energy  level 
fluctuation  for  the  breast  has  a  mean  of  3.36  mm,  which  is 
45%,  larger  than  the  abdominal  wall  correlation  length  of 
2.31  mm.  Consequently,  the  ratio  of  the  average  arrival  time 
fluctuation  effective  correlation  length  to  the  average  energy 
level  fluctuation  correlation  length  is  2.5  for  the  abdominal 
wall  but  only  1.3  for  the  breast  specimens.  This  shows  that 
the  arrival  time  and  energy  level  fluctuation  correlation 
lengths  for  the  breast  are  more  similar  to  each  other  than  are 
those  for  the  abdominal  wall. 

The  focus  obtained  with  uncompensated  breast  tissue 
waveforms  was  predictably  worse  than  that  obtained  with 
abdominal  wall  waveforms  as  is  illustrated  by  the  bar  charts 
in  Fig.  8  for  uncompensated  data.  For  example,  the  mean 
“iO-dB  peripheral  energy  ratio  of  the  focused  breast  data 
was  3.8,  nearly  twice  the  average  value  of  2.1  computed  for 
the  abdominal  data.  Also,  although  not  shown  in  Fig.  8,  the 
mean  —  20-dB  effective  radius  for  the  focused  breast  data 
was  4.6  mm  while  the  corresponding  value  was  3.6  mm  for 
the  abdominal  wall  data.  However,  both  the  average  breast 
and  abdominal  wall  data  sets  first  deviate  10%  from  their 
ideal  effective  radius  curves  at  about  —5  dB. 

Time-shift  compensation  after  backpropagation  im¬ 
proves  the  focus  of  breast  data  more  than  time-shift  compen¬ 
sation  in  the  aperture  does.  However,  in  neither  case  is  the 
average  focus  better  than  that  of  abdominal  wall  data  after 
similar  processing.  The  bar  charts  in  Fig.  8  illustrate  these 
trends.  For  example,  the  average  10%  deviation  level  for 
breast  data  is  reduced  to  —18.4  dB  by  time-shift  compensa¬ 
tion  in  the  aperture  and  to  —22.3  dB  by  time-shift  compen¬ 
sation  after  backpropagation,  but  the  corresponding  values 
for  the  abdominal  wall  data  are  —23.5  and  —26.3  dB,  re¬ 
spectively.  Similarly,  time-shift  compensation  in  the  aperture 
reduces  the  average  -10-dB  peripheral  energy  ratio  to  1.0 


for  the  breast  measurements  and  time-shift  compensation  af¬ 
ter  backpropagation  further  reduces  it  to  0.60,  while  the  cor¬ 
responding  values  for  the  abdominal  wall  measurements  are 
0.63  and  0.47.  Although  both  compensation  techniques  im¬ 
prove  the  focus  of  ultrasonic  waveforms  that  have  propa¬ 
gated  through  breast  tissue,  time-shift  compensation  is  more 
effective  when  performed  after  backpropagation. 

The  distance  of  backpropagation  for  maximum  wave¬ 
form  similarity,  i.e.,  the  position  of  the  equivalent  phase 
screen  in  the  propagation  model  employed  here,  is  indicated 
by  the  data  in  Tables  I  and  II  (as  well  as  in  Tables  III- VI)  to 
be  usually  greater  than  the  specimen  thickness.  However,  as 
already  noted,  a  distance  of  about  8  mm  is  present  between 
the  top  surface  of  the  specimen  and  the  receiving  aperture  in 
the  measurements.  Adding  8  mm  to  the  mean  (±s.d.)  speci¬ 
men  thickness  given  in  Table  II  yields  34.9  (±7.7)  mni  while 
the  mean  backpropagation  distance  is  36.7  (±8.7)  mm.  Thus 
the  data  indicate  that  the  optimum  backpropagation  distance 
is  approximately  the  sum  of  the  specimen  thickness  and  the 
specimen-receiver  separation,  as  in  the  case  of  the  abdominal 
wall  measurements.^^  The  physical  origin  of  this  circum¬ 
stance  is  currently  unknown,  although  observations  and  ex¬ 
periments  noted  earlier  indicate  that  the  origin  is  not  refrac¬ 
tion  at  the  specimen  boundary.  Additional  studies  are  needed 
to  provide  information  about  the  relation  between  tissue  mor¬ 
phology  and  ultrasonic  aberration. 

The  plots  of  arrival  time  and  energy  level  fluctuations  in 
Fig.  5  and  the  statistics  of  Table  III  show  that  reducing  the 
tissue  path  length  decreases  the  severity  of  the  distortion  pro¬ 
duced  in  transmitted  ultrasonic  pulses.  For  both  specimens 
9a  and  11,  the  magnitude  of  the  measured  arrival  time  and 
energy  level  fluctuations  decreases  as  the  tissue  thickness  is 
reduced  while  waveform  similarity  generally  increases.  The 
lack  of  similarity  in  the  patterns  of  fluctuations  as  each  speci¬ 
men  decreases  in  thickness  is  attributed  to  configurational 


TABLE  V.  Breast  tissue  path  statistics  of  wavefront  distortion  for  two  source  positions  12  mm  apart. 


Arrival  time  fluctuations  Energy  level  fluctuations 


Specimen 

number 

Source 

position 

Specimen 

thickness 

(mm) 

rms 

value 

(ns) 

99.5% 

value 

(ns) 

Effective 
corr.  len. 
(mm) 

rms 

value 

(dB) 

99.5% 

value 

(dB) 

Effective 
corr.  len. 
(mm) 

Waveform 

similarity 

factor 

8a 

Left 

30-35 

66.5 

222.0 

4.88 

6.08 

16.35 

4.85 

Right 

30-35 

65.2 

216.0 

4.17 

6.01 

15.11 

4.79 

17 

Left 

20-25 

50.5 

162.3 

2.40 

4.53 

12.46 

2.46 

Right 

20-25 

47.0 

150.4 

2.27 

4.36 

11.93 

2.44 

0.909 
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TABLE  VI.  Breast  tissue  path  focus  characteristics  for  two  source  positions.  The  numbers  in  parentheses  are  values  obtained  at  one  source  position  with 
waveforms  compensated  using  time-shifts  calculated  with  waveforms  from  the  other  source  position,  effective  radius.  PER  ==  peripheral  energy  ratio. 
Unc=uncompensated.  TSC= time-shift  compensation.  BP=backpropagation  followed  by  TSC  BP  Dist^distance  of  backpropagation  for  maximum  waveform 
similarity. 


Specimen 

number 

Source 

position 

BP  dist. 
(mm) 

Unc 

(mm) 

10  dB  r 

TSC 

(mm) 

e 

BP 

(mm) 

Unc 

(mm) 

-20  dB  /• 
TSC 
(mm) 

e 

BP 

(mm) 

Unc 

(mm) 

-30  dB  r 
TSC 
(mm) 

e 

BP 

(mm) 

Unc 

10  dB  PER 

TSC  BP 

Unc 

(dB) 

10%  dev.  lev. 

TSC  BP 

(dB)  (dB) 

8a 

Left 

40 

1.48 

1.09 

(1.13) 

1.06 

(1.12) 

4.39 

1.93 

(2.41) 

1.69 

(2.17) 

4.66 

(6.76) 

3.44 

(6.24) 

2.73 

0.91 

(1.75) 

0.63 

(1.35) 

-1.0 

-12.8 

(-7.7) 

-17.0 

(-8.8) 

Right 

40 

1.43 

1.05 

(1.09) 

1.05 

(1.09) 

4.61 

1.82 

(2.47) 

1.55 

(1.87) 

4.86 

(7.20) 

3.74 

(6.23) 

3.17 

1.05 

(1.91) 

0.67 

(1.34) 

-1.0 

-16.4 

(-11.5) 

-20.6 

(-12.0) 

17 

Left 

30 

1.04 

1.01 

(0.99) 

1.00 

(1.00) 

2.83 

1.43 

(1.76) 

1.40 

(1.67) 

7.52 

4.25 

(6.98) 

2.57 

(5.46) 

3.04 

0.88 

(2.29) 

0.63 

(1.75) 

-13.3 

-24.0 

(-17.7) 

-26.1 

(-19.1) 

Right 

30 

1.06 

1.02 

(1.01) 

1.01 

(1.02) 

2.86 

1.53 

(1.75) 

1.46 

(1.64) 

7.07 

4.67 

(...) 

2.31 

(5.92) 

2.79 

0.88 

(2.32) 

0.63 

(1.78) 

-12.3 

-22.2 

(-16.9) 

-28.1 

(-18.4) 

changes  that  result  from  removal  of  tissue  slabs,  differences 
in  position  of  the  specimen,  and  deformation  from  pressure 
applied  when  the  specimen  was  sliced. 

The  nine  separate  breast  measurements  show  a  similar 
dependence  of  distortion  on  specimen  thickness.  Linear  re¬ 
gression  establishes  significant  relationships  between  breast 
tissue  thickness  and  both  rms  arrival  time  and  energy  level 
fluctuation  values,  with  correlation  coefficients  of  0.729  and 
0.724,  respectively.  The  energy  level  fluctuation  correlation 
lengths  are  somewhat  correlated  to  thickness  as  well.  For  the 
abdominal  wall  specimens,  the  arrival  time  fluctuation  rms 
and  correlation  length  values  are  moderately  correlated  to 
thickness  and  the  waveform  similarity  factor  is  closely  cor¬ 
related,  but  energy  fluctuation  rms  values  and  correlation 
lengths  are  not  correlated  to  thickness.  Since  breast  tissue 
consists  of  a  heterogeneous  arrangement  of  fat,  glandular, 
and  connective  tissue  while  the  abdominal  wall  is  comprised 
of  distinct  fat  and  muscle  layers,  morphological  differences 
are  thought  to  be  the  reason  for  the  observed  differences  in 
fluctuations  produced  by  breast  and  abdominal  wall. 

The  arrival  time  fluctuation  maps  in  Fig.  6  and  focus 
parameters  in  Table  VI  for  measured  wavefronts  with  differ¬ 
ent  source  positions  demonstrate  that  the  time  shifts  that  best 
compensate  the  focus  at  one  position  may  be  less  effective 
when  the  desired  focus  is  moved  12  mm  laterally  and  indi¬ 
cate  that  the  use  of  backpropagation  reduces  this  problem. 
Compensation  in  the  aperture  using  the  time-shifts  calculated 
for  the  adjacent  source  position,  i.e.,  cross  compensation, 
improves  the  focus  in  every  case,  but  is  far  less  effective  than 
compensation  using  the  time-shifts  calculated  for  the  data  set 
being  corrected.  Backpropagation  increases  the  peak  cross¬ 
correlation  coefficient  of  the  arrival  time  fluctuation  map 
pairs  to  0.728  from  0.564  for  specimen  8a  and  to  0.568  from 
0.426  for  specimen  17.  Cross  compensation  is  more  effective 
after  backpropagation,  but  still  does  not  achieve  the  level  of 
focus  improvement  obtained  using  regular  time-shift  com¬ 
pensation  in  the  aperture.  This  implies  a  limitation  on  the 
size  of  the  region  over  which  a  single  pattern  of  compensa¬ 
tion  is  useful.  For  a  related  theoretical  treatment  that  yields 
an  expression  for  the  size  of  the  isoplanatic  patch,  see  Ref. 
19. 
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The  statistics  of  the  arrival  time  distortion  measured  in 
this  study  can  be  compared  to  the  in  vivo  results  of 
Freiburger  et  al  for  a  two-dimensional  array.^  Their  mean 
rms  phase  distortion  value  of  55 .3  ±14.0  ns  obtained  for  35 
scans  on  the  left  breasts  of  seven  volunteers  is  somewhat 
lower  than  the  66.8±12.6  ns  determined  here  for  nine  breast 
specimens.  However,  their  tissue  paths  had  a  mean  length  of 
78.8±21.4  mm,^*^  which  is  much  longer  than  those  measured 
here.  As  discussed  above,  a  longer  path  is  expected  to  pro¬ 
duce  larger  fluctuations  in  arrival  time.  Also,  their  receiver 
elements  were  0.51X3.50  mm^,  which  corresponds  to  a  sur¬ 
face  area  1.7  times  larger  than  that  of  the  elements  used  in 
this  study.  A  larger  measurement  spot  is  associated  with  a 
decrease  in  measured  arrival  time  fluctuations  because  the 
fluctuations  are  averaged  over  the  measurement  spot,  as  de¬ 
tailed  in  a  study^^  that  emulated  the  outputs  of  elements  of 
various  sizes  in  one-dimensional  apertures  from  the  mea¬ 
sured  outputs  of  smaller  elements  in  a  two-dimensional  ap¬ 
erture.  Nevertheless,  the  4.31±1.08-mm  phase  correlation 
length  measured  here  is  similar  to  the  4.21±1.14-mm  value 
obtained  by  Freiburger  et  al^  in  the  azimuthal  direction.  The 
effects  of  the  distortion  measured  in  the  two  cases  on  the 
resulting  focus  are  not  comparable  because  of  substantial  dif¬ 
ferences  in  focusing  parameters. 

The  wavefront  distortion  and  focus  degradation  found  in 
this  study  may  also  be  compared  with  the  already  cited  in 
vivo  results  of  Zhu  and  Steinberg.^"^  Their  observation  of 
severe  amplitude  distortion  accompanying  phase  distortion 
in  narrow-band  measurements  is  analogous  to  energy  level 
fluctuations  accompanying  arrival  time  fluctuations  in  the 
present  wideband  measurements.  The  emphasis  that  Zhu  and 
Steinberg  give  to  the  influence  of  amplitude  distortion  in 
their  analysis  is  strongly  supported  by  the  new  data  given 
here.  Zhu  and  Steinberg  also  showed  representative  recon¬ 
structions  of  source  profiles  that  are  analogous  in  the  present 
study  to  effective  width  curves  of  the  uncompensated  focus 
in  the  array  direction.  The  high  degradation  of  the  source 
profiles  illustrated  in  their  reports  is  qualitatively  similar  to 
the  poor  focus  characteristics  given  here  for  uncompensated 
data.  However,  a  quantitative  comparison  of  source  profile 
and  focus  width  characteristics  is  not  readily  made  with  the 
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FIG.  7.  Comparison  of  breast  and  abdominal  wall  wavefront  distortion  sta¬ 
tistics.  In  each  chart,  the  average  and  standard  deviation  of  the  measure¬ 
ments  within  each  group  are  shown. 


FIG.  8.  Comparison  of  breast  and  abdominal  wall  focus  statistics.  In  each 
chart,  the  average  and  standard  deviation  of  the  measurements  within  each 
group  are  shown.  Uncomp = uncompensated  waveforms.  TSC = time-shift 
compensation.  BP&TSC=backpropagation  before  time-shift  compensation. 
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published  Zhu  and  Steinberg  data  because  they  employed 
different  focus  descriptors,  as  well  as  relatively  narrow-band 
pulses,  a  larger  measurement  spot  size,  and  significantly 
longer  tissue  paths. 

IV.  CONCLUSION 

Ultrasonic  wavefront  distortion  produced  by  transmis¬ 
sion  through  breast  specimens  has  been  investigated  by  re¬ 
cording  and  processing  pulsed  waveforms  to  obtain  maps 
and  statistics  of  arrival  time  and  energy  level  fluctuations. 
The  average  rms  value  of  the  arrival  time  fluctuations  was 
similar  to  that  for  abdominal  wall  specimens  of  like  thick¬ 
ness  while  the  average  rms  value  of  energy  level  fluctuations 
was  greater  than  that  for  similarly  thick  abdominal  wall 
specimens.  The  average  effective  correlation  length  of  the 
arrival  time  fluctuations  for  breast  specimens  was  shorter 
than  that  for  the  abdominal  wall  while  the  average  effective 
correlation  length  of  the  energy  level  fluctuations  was  larger 
than  that  for  the  abdominal  wall.  The  breast  tissue  wave¬ 
forms  showed  more  variation  than  the  abdominal  wall  wave¬ 
forms.  Time-shift  compensation  in  the  aperture  improved  the 
focus  of  pulses  distorted  by  breast  tissue  specimens  and 
time-shift  compensation  after  backpropagation  brought  the 
focus  closer  to  the  ideal,  but  the  improvement  was  generally 
less  than  that  obtained  for  distortion  produced  by  abdominal 
wall.  Repeated  study  of  specimens  with  reduction  of  thick¬ 
ness  showed  that  pulse  arrival  times,  energy  levels,  and  wave 
shape  are  increasingly  altered  as  tissue  path  length  increases. 
Measurements  using  different  source  positions  indicated  that 
time-shift  compensation  in  the  aperture  with  arrival  times 
estimated  when  the  source  was  12  mm  away  firom  the  focal 
position  produced  little  improvement  in  the  focus,  but  that 
backpropagation  made  such  cross  compensation  more  effec¬ 
tive,  although  still  less  effective  than  self-compensation  in 
the  receiving  aperture.  These  results  demonstrate  that  ultra¬ 
sonic  propagation  through  breast  tissue  produces  appreciable 
arrival  time  and  energy  level  fluctuations  and  provide  impor¬ 
tant  new  information  about  the  variety  and  range  of  these 
fluctuations.  They  also  show  that  time-shift  compensation 
improves  the  focus  of  waveforms  distorted  by  breast  tissue, 
especially  when  the  compensation  is  applied  after  wavefront 
backpropagation. 
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ABSTRACT 


Pulse-echo  measurements  from  random  scattering  and  from  a  point  target  have  been 
used  to  quantify  transmitter  beam  size  effects  and  isoplanatic  patch  size  as  well  as  to 
evaluate  the  performance  of  different  aberration  compensation  techniques.  Measurements 
were  made  using  a  single-element  transmitter  with  a  diameter  of  1/2",  1",  or  2",  each  fo¬ 
cused  at  3".  A  tissue-mimicking  scattering  phantom  or  a  point  target  was  used  to  produce 
echoes  that  were  received  in  a  two-dimensional  aperture  synthesized  by  scanning  a  linear 
array.  A  specimen  of  abdominal  wall  was  placed  in  the  reception  path  to  produce  aber¬ 
ration.  B-scan  images  were  formed  with  no  compensation,  with  time-shift  compensation 
in  the  receiving  aperture,  and  with  backpropagation  followed  by  time-shift  compensation. 
The  isoplanatic  patch  size  was  estimated  by  compensating  the  focus  of  a  test  point  target 
with  the  parameters  estimated  for  an  original  point  target  position,  and  observing  the 
deterioration  of  compensation  effects  with  increasing  distance  between  the  test  and  the 
original  point  targets.  The  results  of  the  measurements  using  different  transmitter  di¬ 
ameters  quantify  the  improvement  of  time-delay  estimation  with  the  increase  in  wavefront 
coherence  that  accompanies  decreased  transmitter  beam  size.  For  7  specimens,  the  average 
isoplanatic  patch  size  determined  from  a  10%  increase  in  the  -10  dB  effective  diameter 
was  16.7  mm  in  the  azimuthal  direction  and  39.0  mm  in  the  range  direction.  These  sizes 
increased  after  backpropagation  to  19.0  mm  and  41.4  mm,  respectively.  For  the  1/2",  1", 
and  2"  diameter  transmitters,  the  average  contrast  ratio  improvement  was  2.0  dB,  2.1  dB, 
and  2.8  dB,  respectively,  with  time-shift  compensation,  and  2.3  dB,  2.7  dB,  and  3.5  dB, 
respectively,  with  backpropagation  of  20  mm  followed  by  time-delay  estimation  and  com¬ 
pensation.  The  investigation  indicates  that  a  tightly  focused  transmitter  beam  is  necessary 
to  create  a  scattered  wavefront  satisfactory  for  time-shift  estimation,  the  isoplanatic  patch 
is  about  twice  as  long  in  the  range  direction  as  in  the  azimuthal  direction,  and  backprop¬ 
agation  followed  by  time-shift  compensation  provides  better  compensation  of  distortion 
than  time-shift  compensation  alone. 
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INTRODUCTION 


The  estimation  and  correction  of  ultrasonic  wavefront  distortion  caused  by  propaga¬ 
tion  through  an  inhomogeneous  medium  (as  most  tissues  are)  has  been  studied  by  academic 
and  industrial  researchers  using  various  approaches.  Some  started  from  data  measured  with 
one-dimensional,  clinical  arrays  having  a  large  elevation.  Some  developed  algorithms  and 
systems  with  emphasis  on  real-time  implementation.  Some  measured  the  distortion  effects 
using  a  relatively  controlled  yet  clinically  relevant  configmation.  This  latter  approach  has 
yielded  transmission  measurements  using  a  spherical  wave  source  and  two-dimensional  re¬ 
ceiving  apertme  with  a  small  (about  1  mm^)  measurement  spot  size  and  has  provided 
basic  data  that  describe  the  characteristics  of  wavefront  aberration  caused  by  tissue.  The 
next  logical  step  is  to  apply  this  information  to  correct  for  tissue  aberration  in  a  pulse-echo 
setting  using  random  scattering,  similar  to  conditions  commonly  encountered  in  diagnostic 
imaging.  A  study  that  does  so  and  aims  to  provide  considerable  insight  into  aberration 
correction  in  the  pulse-echo  mode  is  reported  here. 

Although  the  existence  of  wavefront  distortion  has  been  known  since  almost  the  be¬ 
ginning  of  b-scan  imaging  [1,2],  the  need  for  compensation  of  this  distortion  has  only 
recently  become  great  with  the  advent  of  large,  high-frequency  arrays.  Extensive  investi¬ 
gation  of  aberration  measurement  and  correction  conducted  in  view  of  this  need  started 
about  10  years  ago.  Flax  and  O’Donnell  [3]  measured  the  time  distortion  based  on  the 
cross-correlation  of  waveforms  received  by  neighboring  elements  of  a  linear  array.  Waag,  .e^ 
al.  [4,5]  reported  measurements  and  analysis  of  wavefront  distortion  caused  by  propagation 
through  abdominal  wall  and  calf  liver.  Nock,  et  al.  [6]  proposed  a  maximum  brightness 
method  that  was  implemented  in  real-time.  Focus  compensation  in  these  early  studies 
concentrated  on  the  correction  of  arrival  time  fluctuations.  However,  later  studies  [7-11] 
have  shown  that  correction  of  arrival  time  fluctuation  alone  is  not  sufficient  to  compen¬ 
sate  fully  for  ultrasonic  distortion  caused  by  propagation  through  an  extended  thickness 
of  tissue.  Measurements  [7,8]  indicate  that  the  propagated  wavefront  contains  not  only 
time-shift  aberration  but  also  substantial  energy  level  fluctuation  and  significant  waveform 
changes.  Focusing  studies  [9,10]  show  that  time-shift  compensation  improves  mainly  the 
point  resolution  but  that  the  sidelobes  remain  considerably  higher  than  that  of  an  ideal  fo¬ 
cus.  Based  on  these  results  and  observations,  an  improved  wavefront  distortion  model  that 
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consists  of  a  time-shifting  screen  placed  some  unknown  distance  away  from  the  measure¬ 
ment  aperture  was  proposed  and  employed  with  measured  data  [12],  The  results  obtained 
using  this  model  with  abdominal-wall  data  [12]  and  breast  data  [9]  suggest  that  a  better 
way  to  compensate  for  wavefront  distortion  is  first  to  backpropagate  the  wavefront  by  an 
appropriate  distance,  and  then  to  time-shift  compensate  the  wavefront  at  that  distance. 

Accurate  estimation  of  the  distortion  is  required  to  model  and  compensate  the  dis¬ 
tortion  effectively.  This  estimation,  which  includes  calculation  of  arrival  time  fiuctuation 
from  wavefront  data,  is  a  difficult  task,  since  severe  waveform  distortion  can  invalidate 
the  usual  definition  and  estimation  of  arrival  time.  Additionally,  for  randomly  distributed 
scatterers  illuminated  by  a  focused  beam,  the  scattered  signals  are  known  to  decorrelate 
with  increasing  spatial  distance  between  the  measurement  spots.  The  degree  of  decorrela¬ 
tion  is  related  to  the  illumination  beam  size,  and  is  governed  by  the  van  Cittert-Zernike 
theorem  [13,14]  under  the  assumptions  of  incoherent  scattering  and  a  homogeneous  prop¬ 
agation  path.  Another  factor  that  influences  the  perceived  distortion  effect  is  the  size 
of  receiver  elements.  Simulations  of  a  one-dimensional  aperture  using  measurements  in 
a  two-dimensional  aperture  have  shown  [15]  that  accurate  measurement  of  distortion  re¬ 
quires  small  size  elements  to  avoid  significant  spatial  averaging  over  the  element  surface. 
All  the  foregoing  points  must  be  addressed  for  the  appropriate  measurement,  estimation, 
and  correction  of  wavefront  distortion. 

The  application  of  wavefront  aberration  correction  in  in  vivo  clinical  diagnostic  imag¬ 
ing  requires  the  estimation  of  a  distortion  profile  using  randomly  scattered  waveforms 
produced  by  a  beam  of  illumination.  The  distortion  effect  of  the  inhomogeneous  medium 
may  vary  for  different  regions  of  the  scattering  medium  because  wavefronts  originating 
from  different  positions  propagate  through  different  paths  in  the  inhomogeneous  medium 
and  are  distorted  differently.  Critical  questions  for  the  practical  implementation  of  wave- 
front  aberration  correction  are  the  following.  How  accurately  can  the  distortion  profile 
be  estimated?  For  how  large  a  spatial  region  is  a  single  estimation  suitable?  How  much 
improvement  can  be  expected  by  distortion  compensation?  The  answer  to  these  questions 
naturally  depends  on  the  severity  of  the  distortion  and  the  compensation  technique  being 
used.  In  this  paper,  pulse-echo  experiments  and  calculations  designed  to  address  the  above 
questions  are  described  to  provide  quantitative  results  about  the  following. 
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1.  The  effect  of  the  transmitter  beam  size  on  the  coherence  of  scattered  wavefront  and 
on  the  accuracy  of  time-delay  estimation. 

2.  The  isoplanatic  patch  size,  ie.,  the  size  of  a  region  for  which  the  distortion  effect  can 
be  corrected  using  the  same  set  of  parameters. 

3.  The  use  of  backpropagation  in  the  pulse-echo  mode  using  randomly  scattered  wave¬ 
forms,  and  the  effectiveness  of  time-shift  compensation  with  and  without  backpropa¬ 
gation. 

The  organization  of  this  paper  is  the  following.  First,  data  acquisition  is  described 
along  with  the  considerations  of  experimental  design.  Next,  issues  in  data  processing, 
particularly  in  time-delay  estimation,  are  discussed,  and  the  formation  of  b-scan  images 
with  dynamic  receive  focusing  is  described  for  the  particular  geometry  used  in  the  mea¬ 
surements.  Also,  methods  used  for  the  evaluation  of  wavefront  coherence,  image  contrast 
ratio,  and  point  focus  are  described.  Then,  results  from  seven  specimens  of  abdominal 
wall  are  presented  for  waveform  similarity  factor,  rms  arrival  time  fluctuation,  isoplanatic 
patch  size,  and  contrast  ratio.  The  results  are  next  discussed  and  Anally  conclusions  are 
drawn. 
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I.  DATA  ACQUISITION 


A.  Measurement  Configuration 

The  measurement  apparatus  consisted  of  a  transmitter,  a  scattering  phantom,  a  speci¬ 
men,  and  a  receiving  linear  array,  configured  as  shown  in  Fig.  1.  The  apparatus  is  immersed 
in  distilled  water  with  temperatme  maintained  at  30°  C  using  a  heater  and  circulation.  In 
this  configuration,  the  transmit  beam  is  unperturbed  and  only  the  receive  beam  is  dis¬ 
torted  by  the  aberration  medium,  for  which  specimens  of  human  abdominal  wall  were 
used.  The  arrangement  permits  control  of  the  transmit  beam  and  compensation  of  the 
receive  beam  by  different  techniques  using  data  measured  in  a  two-dimensional  aperture. 
Two  scan  movements  used  in  the  measurement  are  illustrated  in  Fig.  2.  The  first  movement 
scans  the  linear  array  to  form  a  two-dimensional  receiving  aperture.  The  second  move¬ 
ment,  analogous  to  conventional  b-scan  imaging,  scans  the  transmitter  and  the  receiving 
aperture  together  to  form  different  scan  lines. 

The  transmit  beam  was  produced  using  a  single  element  transducer  and  the  beam 
size  was  controlled  by  using  transducers  of  different  diameters,  i.e.,  1/2",  1",  and  2",  each 
focused  at  3".  The  fixture  was  designed  so  that  the  focal  point  of  each  transducer  was 
at  the  same  spatial  point.  Due  to  the  different  electrical  impedances  of  the  transducers, 
a  custom-made  pulser  with  individual  impedance  matching  networks  for  each  transducer 
was  used  to  msiximize  the  energy  radiated  into  the  mediiun. 

A  tissue-mimicking  phantom  with  a  10  mm-diameter  scatterer-free  cylinder  to  simulate 
a  blood  vessel  or  a  cyst  was  used  as  the  scattering  medium.  To  increase  scattering  signal 
strength  and  reduce  attenuation  in  the  phantom,  the  phantom  was  made  of  agar  (40  grams 
per  liter  of  water)  and  glass  spheres  of  diameters  around  75~90  jum  (4  grams  per  liter  of 
agar  solution).  This  resulted  in  about  10  dB  greater  signal  strength  compared  to  previously 
used  graphite-gel  phantoms.  The  final  signal-to-noise  ratio  achieved  with  the  measurement 
system  averaged  21.9  dB,  22.9  dB,  and  19.0  dB  for  seven  specimens  of  abdominal  wall  using 
the  1/2",  1",  and  2"  transducer,  respectively. 

In  place  of  the  phantom,  a  point  target  was  positioned  at  the  focal  point  of  the  trans¬ 
mitting  transducer  to  produce  a  spherical  scattered  wavefront  with  the  highest  possible 
spatial  coherence  so  that  time  delays  estimated  from  such  a  wavefront  could  be  used  as  a 
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standard  against  which  time  delays  estimated  using  random  scattering  from  the  phantom 
could  be  compared. 

A  diagram  of  the  data  acquisition  system  is  shown  in  Fig.  3.  A  free-running  oscillator, 
synchronized  with  the  20  MHz  clock  signal  of  the  analog-to-digital  (A/D)  converter,  trig¬ 
gered  a  custom-made  pulser  that  excited  the  single-element  transducers  through  impedance 
matching  networks  individually  tuned  for  each  transducer.  The  scattered  signal  was  re¬ 
ceived  by  a  128-element  linear  array  and  the  elements  were  individually  accessed  using 
a  multiplexer.  The  pitch  of  the  array  is  0.72  mm  and  the  element  size  in  the  elevation 
direction  was  limited  to  1.00  mm  to  avoid  excessive  spatial  averaging  across  the  element 
surface.  The  output  from  the  accessed  element  then  passed  through  a  TGC  amplifier  and 
a  fixed  gain  amplifier.  The  signal  was  next  A/D  converted  with  a  precision  of  12  bits  at  a 
rate  of  20  MHz.  The  measurement  at  each  element  was  repeated  6  times  (determined  by 
hardware)  so  that  data  averaging  could  be  used  to  reduce  electronic  noise  and  the  results 
were  stored  in  a  waveform  memory.  After  data  were  acquired  for  all  elements  at  each 
elevation  and  scan  position,  the  data  were  transferred  to  a  PC  via  GPIB,  averaged,  and 
stored  on  a  disk  temporarily  before  being  sent  (via  Internet)  for  subsequent  processing  on 
an  IBM  SP2  computer.  Parameters  of  the  measurement  configuration  are  tabulated  in 
Table  I.  As  indicated  in  Table  I  and  shown  in  Fig.  2,  the  linear  array  was  moved  over  12 
elevation  positions  to  form  the  two-dimensional  receiving  aperture,  and  the  transmitter 
and  linear  array  assembly  was  moved  over  20  positions  to  obtain  20  scan  lines.  Because 
the  scan  movement  was  mechanical,  the  scan  step  size  could  be  chosen  independent  of  the 
element  pitch  of  the  linear  array.  The  size  of  each  dataset  was  (2  bytes/sample)  x  (776 
samples/element)  x  (66  elements/elevation)  x  (12  elevations/scan)  x  (20  scans),  or  about 
25  MBytes. 

A  water-path  point  target  echo  waveform  and  its  amplitude  spectrum  are  presented 
in  Fig.  4  to  show  the  center  frequency  and  bandwidth  of  the  measurement  system.  Typi¬ 
cal  pulse-echo  wavefronts  obtained  using  a  point  target  and  using  the  random  scattering 
phantom  are  shown  in  Fig.  5  for  a  representative  tissue  path. 

B.  Measurements  for  Evaluating  Isoplanatic  Patch  Size 

The  point  target  positioned  at  the  center  of  the  transmit  focus  provided  a  means  for 
measuring  the  isoplanatic  patch  size,  ie.,  the  size  of  the  scattering  region  for  which  the 
distortion  effect  of  the  inhomogeneous  medium  is  about  the  same. 
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Data  for  the  evaluation  of  isoplanatic  patch  size  in  the  range  direction  were  obtained 
by  moving  the  transmitter/point-target  assembly  in  the  vertical  direction  manually  and 
measuring  the  scattered  wavefront  in  a  two-dimensional  aperture  for  each  point  target 
position.  See  Fig.  6.  The  topmost  point  target  position  was  about  81  mm  below  the 
apertme,  and  5  sequential  positions  were  measured  in  an  increment  of  7  mm,  yielding  5 
sets  of  wavefront  data. 

The  sets  of  wavefront  data  were  independently  corrected  for  geometric  path  length 
differences,  after  which  all  the  wavefronts  were  essentially  flat  having  only  small  fluctua¬ 
tions  from  wavefront  distortion.  The  geometric  correction  was  performed  by  estimating  an 
arrival  time  surface,  fitting  a  spherical  surface  to  the  arrival  time  surface,  and  shifting  the 
waveforms  with  the  fitted  delays.  Next,  a  time  delay  map  was  estimated  for  the  wavefront 
corresponding  to  the  topmost  point  target,  and  this  was  used  to  time  shift  each  of  the 
5  wavefronts.  In  this  way,  the  wavefront  corresponding  to  the  topmost  point  target  was 
time-shift  compensated  with  delays  estimated  from  itself,  which  is  called  here  self  compen¬ 
sation,  and  the  other  wavefronts  were  time-shift  compensated  with  delays  estimated  for  a 
different  point  target  position,  which  is  called  here  cross  compensation.  In  cross  compen¬ 
sation,  a  test  point  target  wavefront  is  said  for  brevity  to  be  compensated  with  distortion 
parameters  estimated  for  the  original  point  target.  The  compensated  wavefronts  were  then 
each  focused  to  produce  an  image  of  the  point  target. 

The  data  for  the  evaluation  of  isoplanatic  patch  size  in  the  azimuthal  direction  was 
similarly  obtained  by  scanning  the  transmitter  and  point-target  assembly  in  the  horizontal 
or  array  direction.  During  the  scan,  the  specimen  was  stationary  while  the  receiving  array 
translated  with  the  transmitter.  This  resulted  in  a  relative  shift  between  the  specimen  and 
the  two-dimensional  receiving  aperture.  In  order  to  maintain  interrogation  of  the  same 
tissue  path,  a  different  portion  of  the  array  was  used  for  forming  the  receiving  aperture. 
Therefore,  the  azimuthal  scan  step  size  was  made  to  be  the  same  as  the  element  pitch 
in  the  linear  array,  i.e.,  0.72  mm,  and  with  each  scan,  the  portion  of  array  used  to  form 
the  receiving  aperture  was  shifted  by  one  element  in  the  direction  opposite  to  the  scan 
movement.  The  number  of  such  scans  was  20.  The  vertical  distance  from  the  point  target 
to  the  aperture  in  these  measurements  was  about  95  mm.  In  processing,  the  20  measured 
wavefronts  were  each  corrected  for  geometric  path  length  difference,  and  a  time-shift  map 


8 


was  estimated  for  the  first  wavefront.  This  map  was  used  to  time-shift  compensate  each 
of  the  20  wavefronts.  The  time-shift  compensated  wavefronts  were  then  each  focused  to 
produce  an  image  of  the  point  target. 
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II.  DATA  PROCESSING 


A.  Time-Delay  Estimation 

Time-delay  estimation  is  the  central  issue  in  data  processing.  If  time  delays  are  not 
estimated  correctly,  the  full  potential  improvement  that  results  from  time-shift  compen¬ 
sation  (with  or  without  backpropagation)  cannot  be  achieved,  and  the  effectiveness  of 
techniques  that  employ  time-shift  compensation  cannot  be  properly  assessed.  Similarly, 
for  the  measurement  of  isoplanatic  patch  size,  if  the  point  target  focus  is  not  optimal  af¬ 
ter  compensation  with  time  delays  estimated  for  that  position,  then,  as  the  point  target 
moves  away,  the  deterioration  of  the  compensation  effect  will  be  less  noticeable,  leading  to 
an  overestimation  of  the  isoplanatic  patch  size. 


The  problem  of  time-delay  estimation  for  a  two-dimensional  array  of  MxN  waveforms 
Si{t)  for  i  =  1  MN  can  be  posed  as  an  optimization  problem  in  which  the  maximization 


of 
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is  sought  subject  to  a  smoothness  requirement  on  the  map  of  Tj’s.  The  quantity  J  can  be 
interpreted  as  the  energy  of  the  coherent  sum  of  the  aligned  waveforms,  and  is  proportional 
to  the  waveform  similarity  factor  defined  in  Ref.  [12]  and  the  focusing  criterion  defined 
in  Ref.  [16].  The  smoothness  requirement  on  the  r*  map  comes  from  the  observation 
that,  since  the  time  delay  aberration  is  induced  by  a  specimen  of  soft  tissue  with  smooth 
structures,  a  two-dimensional  map  of  arrival  time  is  not  likely  to  contain  any  abrupt  jumps 
or  local  points  that  are  very  different  from  their  neighboring  points.  This  smoothness 
requirement  is  experimentally  supported  by  the  smoothness  of  point  target  wavefronts 
measured  through  tissue  specimens.  Although  direct  optimization  of  J  is  computationally 
intensive  and  difficult  because  numerous  local  minima  exist,  this  general  point  of  view  was 
used  to  guide  the  development  of  specific  time  delay  estimation  algorithms  for  the  study 
reported  here. 


In  the  current  study,  time  delays  were  estimated  from  two  types  of  wavefront:  point 
target  wavefronts  and  random  scattering  wavefronts.  These  two  types  of  wavefront  have 
different  spatial  coherence  properties  in  addition  to  different  temporal  waveform  shapes. 
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The  point  target  wavefront  consists  of  a  single  pulse  and  is  coherent  throughout  the  aper¬ 
ture  (that  is,  in  the  absence  of  propagation  aberration,  all  the  pulses  are  the  same),  whereas 
the  scattering  wavefront  consists  of  random  waveforms  that  lasts  a  long  time  and  has  a 
spatial  coherence  that  decreases  with  increasing  separation  between  points  in  the  aperture. 
Previously,  the  authors  have  proposed  two  different  algorithms  for  time-delay  estimation 
in  a  two-dimensional  aperture.  These  are  the  least-mean-square-error  method  [10]  and 
the  reference  waveform  method  [12].  Both  methods  employed  cross-correlation  followed 
by  spline  interpolation  to  detect  the  peak  as  the  basic  means  to  find  the  relative  time 
delay  between  two  given  waveforms.  However,  the  reference  waveform  method  assumes 
ail  the  waveforms  in  the  aperture  to  be  the  same,  and  so  is  only  suitable  for  point  target 
wavefronts.  Another  method  known  in  the  literature  for  time-delay  estimation  in  a  two- 
dimensional  aperture  is  the  row-sum  method  [17],  but  this  method  in  its  presented  form 
also  assumes  the  waveforms  to  be  similar  throughout  each  row  and  is  considered  to  be 
suitable  only  for  point  target  wavefronts. 

The  least-mean-square-error  method  uses  relative  time  delays  for  neighboring  wave¬ 
forms  and  requires  the  waveforms  be  only  locally  coherent.  This  method  is  suited  for 
estimating  the  arrival  time  surface  for  a  wavefront  that  arises  from  random  scattering. 
The  difficulty  with  this  method  is  the  need  to  determine  an  arrival  time  when  neighboring 
waveforms  are  poorly  correlated.  In  the  current  study,  the  least-mean-square-error  method 
was  improved  by  using  reference  points  to  guide  the  search  for  peak  positions  in  the  cross¬ 
correlation  function  and  using  weights  in  the  solution  of  the  least-squares  problem.  The 
improved  version  was  used  as  the  method  for  time-delay  estimation. 

The  least-mean-square-error  method  for  time  delay  estimation  in  a  two-dimensional 
aperture  has  been  described  elsewhere  [10],  so  only  an  outline  is  provided  here.  In  this 
method,  for  a  rectangular  array  of  waveforms,  relative  time  delays  are  found  for  neighboring 
waveforms  in  the  horizontal,  vertical,  and  the  two  diagonal  directions.  These  relative  time 
delays  are  then  used  to  write  a  set  of  linear  equations  of  the  form  Tj  —  r,-  =  dij,  where 
Tj  represents  the  unknown  time  delay  for  waveform  i  and  dij  represents  the  given  relative 
time  delay  between  waveforms  i  and  j.  For  an  array  oi  MxN  waveforms,  the  number  of 
such  equations  is  4MN  —  3(M  +  N)  +  2,  and  the  number  of  unknowns  is  MN  —  1.  (One 
waveform  is  used  as  a  time  reference.)  Therefore,  the  number  of  equations  is  approximately 
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four  times  the  number  of  unknowns.  These  equations  are  solved  in  the  least-squares  sense 
by  minimizing 

-  Tj  -  dijf]  , 

»  3 

where  the  summations  for  i  and  j  are  over  all  the  relative  time  delays  dij's  that  have  been 
given,  and  Wij  is  the  weight  for  dij.  The  solution  of  this  least-squares  problem  can  exploit 
the  simple  structure  of  the  equations  and  is  computationally  inexpensive  compared  to  the 
calculation  of  all  the  cross-correlation  functions  that  are  needed  to  find  the  relative  time 
delays. 

The  estimation  of  relative  time  delay  between  two  waveforms  was  b£ised  on  the  position 
of  a  peak  in  the  cross-correlation  function  of  the  two  waveforms.  However,  instead  of 
seeking  the  peak  position  blindly,  a  reference  point  was  used  to  guide  the  search  for  the 
peak  position  and  to  provide  a  weight  for  the  obtained  result.  To  reduce  computation, 
the  cross-correlation  function  of  neighboring  waveforms  was  calculated  over  only  ±10  lags. 
The  search  in  the  cross-correlation  function  returned  a  peak  position  that  was  closest  to 
the  reference  point  and  that  position  was  used  for  the  relative  time  delay. 

The  weight  of  the  relative  time  delay  so  obtained  was  calculated  using  the  distances 
of  the  reference  point  to  the  two  closest  peaks  on  either  side  of  the  reference  point.  When 
the  two  distances  were  about  equal,  meaning  which  peak  to  choose  was  not  clear,  a  low 
confidence  was  placed  on  the  result.  When  one  peak  was  much  closer  than  the  other  to 
the  reference  point,  a  high  confidence  was  placed  on  the  result.  The  particular  weight 
function  used  is  log[(do  ±  e)/{db  ±  e)],  where  da  is  the  larger  and  dj,  is  the  smaller  of  the 
two  distances  from  the  reference  point  to  the  two  nearest  peaks,  and  e  is  a  small  positive 
number  added  to  avoid  singularity  in  the  calculation. 

The  reference  value  initially  used  to  guide  the  peak  search  was  zero  for  every  pair, 
meaning  that  in  the  absence  of  any  a  priori  knowledge,  the  time  difference  of  neighboring 
waveforms  was  assumed  to  be  zero.  This  is  a  reasonable  initial  guess  so  long  as  the 
wavefront  has  already  been  compensated  for  geometric  path  length  differences.  After  the 
least-squares  problem  was  solved,  a  new  set  of  references  was  obtained  using  the  current 
solution,  i.e.,  r*  —  tj  was  used  as  the  reference  value  for  finding  a  new  set  of  dij.  With 
the  newly  obtained  dij,  the  least-squares  problem  was  solved  again  to  yield  new  values  for 
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Ti’s.  This  was  repeated  until  no  dij  was  altered,  and  the  final  Tj’s  were  used  as  the  result 
of  time-delay  estimation. 


B.  Dynamic  Receive  Focusing 

Dynamic  focusing  in  reception  was  performed  along  the  axis  of  the  transmit  beam. 
The  basic  idea,  like  that  of  dynamic  focusing  in  present  clinical  systems,  is  to  have  the 
receive  focus  track  the  transmit  wavefront.  The  complication  in  this  study  is  that  the 
transmit  and  receive  apertures  did  not  overlap  and  an  angle  of  about  35°  existed  between 
the  transmit  axis  and  the  normal  of  the  receiving  aperture.  For  tracking  the  transmit 
wavefront  in  receive  focusing,  the  position  of  the  transmit  wavefront  at  each  instant  was 
estimated  by  fitting  a  spherical  surface  to  the  estimated  arrival  time  surface,  assuming  that 
the  sound  speed  in  the  scattering  medium  is  known.  A  coordinate  system  was  established 
on  the  surface  of  the  two-dimensional  receive  aperture  in  which  x-axis  was  along  the 
array  direction,  y-axis  was  along  the  elevation  direction,  and  z-axis  was  normal  to  the 
aperture  surface,  as  diagramed  in  Fig.  7.  The  fitting  algorithm  used  to  find  the  center 
of  a  given  spherical  arrival  time  surface  r(x,  y)  is  described  in  an  appendix.  For  random 
scattering,  the  arrival  time  was  estimated  using  the  portion  of  echoes  that  arose  from 
the  focal  region  of  the  transmission  and  had  the  highest  spatial  coherence.  The  fitted 
parameters  (Xq,  Yq,  Zq)  then  represented  the  spatial  origin  for  that  portion  of  echoes  and 
corresponded  to  the  position  of  the  transmit  focal  point.  To  calculate  the  spatial  origin 
for  other  portions  of  echoes,  the  angle  a  of  the  transmit  axis  with  respect  to  the  normal 
of  the  receiving  aperture  was  used.  As  illustrated  in  Fig.  7,  a  distance  I  (At)  along  the 
transmission  axis  was  determined  from  the  focal  point  such  that  the  propagation  time 
along  the  new  path  is  At  longer  than  that  for  the  path  that  goes  through  the  focal  point. 
The  spatial  origin  of  the  wavefront  that  arrives  At  later  than  the  focal  point  wavefront  is 
then  given  by 


(xo(t),yo(t),zo(t))  =  (Xo,Vo  +  l(At)sina,Zo  +  l(At)coso!) . 


To  bring  a  wavefront  originating  from  (xo(t),yo(t),zo(t))  into  focus,  the  waveform  Si(t) 
received  at  (xi,yi,0)  was  shifted  in  time  by  the  amount 


Vi(t)  =■- 


^ (xo(t)  -  Xi)2  -t-  (yo(t)  -  Pi)^  -h  z^(t)  -  \l xl(t) yl(t)  ^  zl(t) 
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This  transformed  the  waveform  Si(t)  into  s^(t)  given  by 

Since  the  time-shift  rii{t)  is  not  a  constant  but  a  function  of  time,  the  waveforms  were  not 
only  time-shifted  but  also  stretched  or  compressed  after  dynamic  focusing.  To  avoid  unde¬ 
sirable  time  quantization  effects,  spline  interpolation  was  used  to  resample  the  waveform 
so  that  need  not  be  quantized  to  be  an  integral  number  of  sampling  intervals. 

C.  Single  Transmit  Image  Formation 

Two  kinds  of  images  were  formed  with  the  data.  The  first  kind  is  similar  to  the  so- 
called  explososcan  imaging  [18]  in  which  a  single  transmit  beam  illuminates  the  medium 
and  the  received  data  are  processed  to  yield  an  image  of  the  medium  with  the  given 
illumination.  The  second  kind  mimics  regular  b-scan  imaging  in  which  both  the  transmit 
and  receive  beams  are  scanned  simultaneously  and  a  line  in  the  b-scan  is  obtained  for  each 
scan  position. 

The  first  step  common  to  both  kinds  of  image  formation  was  the  estimation  of  a 
time-delay  surface  using  the  portion  of  wavefront  that  arises  from  the  transmit  focal  point 
and  has  the  highest  spatial  coherence.  The  signal  length  used  for  time-delay  estimation  is 
an  important  parameter.  Since  the  signals  arise  from  random  scattering,  the  signals  are 
themselves  random.  Use  of  a  longer  signal  length  helps  overcome  the  effects  of  randomness 
on  the  estimated  time  delays.  However,  a  longer  signal  extends  the  region  in  the  scattering 
medium  from  which  the  wavefronts  originate  and  the  effective  distortion  profile  of  the 
specimen  may  consequently  change.  Also,  the  spatial  coherence  property  of  the  wavefronts 
scattered  from  different  positions  may  change  because  the  transmit  beam  has  a  finite  depth 
of  field.  Therefore,  a  compromise  was  used  to  determine  the  signal  length  for  time-delay 
estimation.  A  signal  length  of  128  samples,  corresponding  to  a  segment  length  of  about 
5  mm,  was  chosen  and  found  to  be  appropriate. 

The  time-delay  surface  was  estimated  by  first  removing  an  approximate  geometric 
arrival  time  surface  from  the  data  and  then  applying  the  least-mean-square-error  method 
to  the  shifted  wavefront.  The  approximate  shifts  were  then  added  to  the  least-mean-square- 
error  solution  to  give  the  total  time-delay  surface.  The  reason  for  this  two-stage  process  is 
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that  the  current  implementation  of  the  least-mean-square-error  method  expects  a  wavefront 
that  does  not  contain  significant  curvature.  Next,  a  spherical  fit  to  the  time-delay  surface 
was  used  to  find  the  spatial  position  of  the  focal  point  relative  to  the  receiving  aperture. 
The  fitted  spherical  centers  appeared  to  coincide  with  the  expected  transmit  focal  point 
positions.  Then,  dynamic  receive  focusing  shifts  were  applied  to  the  received  wavefiront 
using  information  about  the  spatial  origin  of  the  echoes  obtained  as  described  previously. 
The  resulting  wavefront  no  longer  contained  any  spherical  curvature,  and  was  essentially 
flat  in  the  absence  of  aberration.  This  intermediate  wavefront,  called  here  the  dynamically 
shifted  wavefront,  plays  an  important  role  in  subsequent  processing. 

A  single  transmit  image,  ie.,  an  image  formed  from  a  single  transmit  beam,  is  similar 
to  a  snap-shot  picture  of  the  medium.  In  the  processing,  the  dynamically  shifted  wavefront 
was  windowed  spatially  using  a  Hamming  window,  passed  through  an  imaging  lens  of  fixed 
focal  length,  and  subsequently  propagated  to  the  focal  plane.  The  focal  plane  time  history 
then  yielded  an  image  of  the  medium  under  the  given  illumination.  This  image  can  also 
be  interpreted  as  the  convolution  of  the  transmit  beam  with  the  receive  beam  because 
the  image  is  obtained  by  fixing  the  transmit  beam  and  steering  the  receive  beam.  The 
propagation  from  the  apertme  plane  to  the  focal  plane  was  performed  by  the  shift-and-add 
method  using  sphne  interpolation  to  implement  arbitrary  time-shifts  of  the  waveforms.  For 
the  focal  length  of  the  imaging  lens,  the  distance  of  the  transmit  focal  point  to  the  origin 
in  the  receiving  aperture  was  used.  Although  the  focal  plane  time  history  corresponds  to  a 
3-D  volume,  which  implies  that  a  volume  image  of  the  medium  can  be  obtained,  attention 
in  the  current  study  was  limited  to  the  time-history  along  the  central  elevation  so  the 
result  is  an  image  of  a  single  cross  section. 

D.  B-Scan  Image  Formation 

B-scan  images  were  formed  without  compensation  (UNC),  with  time-shift  compensa¬ 
tion  (TSC),  and  with  backpropagation  followed  by  time-shift  compensation  (BP-hTSC). 

To  form  an  uncompensated  b-scan  image,  the  waveforms  in  the  dynamically  shifted 
wavefront  were  simply  added  together.  This  represents  a  focusing  operation  by  an  imaging 
lens  and  the  result  corresponds  to  the  time  history  at  the  focal  point  of  the  imaging  lens. 
The  Hilbert  transform  [19]  was  used  to  find  the  analytic  envelope  of  the  waveform  resulting 
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from  the  summation.  The  analytic  envelopes  from  sequential  scan  lines  were  arranged  side- 
by-side  to  form  a  b-scan  image. 

To  form  a  time-shift  compensated  b-scan  image,  a  time-delay  map  was  estimated  for 
each  dynamically  shifted  wavefront,  and  the  result  was  used  to  time-shift  the  wavefront 
at  each  scan  position  individually.  After  such  time-shifting,  the  waveforms  in  the  wave- 
front  are  better  aligned  in  time,  and  summation  results  in  a  signal  of  larger  amplitude. 
For  wavefronts  obtained  with  the  1/2"  transmitter,  time-shift  compensation  was  also  per¬ 
formed  using  the  time-delays  estimated  from  the  corresponding  2"  data.  The  subsequent 
processing  was  the  same  as  that  for  the  uncompensated  case. 

To  form  an  image  using  backpropagation  followed  by  time-shift  compensation,  the 
dynamically  shifted  wavefront  was  first  backpropagated  by  a  fixed  distance  of  20  mm.  The 
distance  was  fixed  because  of  the  large  amount  of  computation  needed  to  carry  out  the 
backpropagation  and  also  because  of  the  previous  finding  [12]  that  the  result  is  relatively 
insensitive  to  the  exact  distance  of  backpropagation.  In  addition,  trial  computations  con¬ 
firmed  that  the  compensation  effect  did  not  change  much  with  the  present  data  when  the 
backpropagation  distance  was  varied  by  a  few  millimeters.  Since  the  specimen  thickness 
ranged  from  10  to  35  mm  and  a  distance  of  5  to  8  mm  existed  between  the  receiving  aper¬ 
ture  and  the  top  surface  of  the  specimen,  backpropagation  by  20  mm  brought  the  wavefront 
about  half-way  into  the  specimen.  The  backpropagation  was  implemented  using  the  angu¬ 
lar  spectrum  method  [20].  In  the  calculation,  an  FFT  was  used  to  calculate  the  temporal 
as  well  as  the  spatial  Fourier  transform.  Waveforms  of  776  samples  were  zero-padded  to 
1024  samples,  and  aperture  size  of  12  x  66  elements  was  zero-padded  to  16  x  128.  Since 
the  same  backpropagation  function  of  the  form  exp  {J2'Kz^Jl/X^  -  fx  ~  fy^  was  needed 
to  backpropagate  the  wavefront  for  each  scan  line,  the  backpropagation  function  was  com¬ 
puted  for  only  one  scan  line  and  then  stored  for  use  with  other  scan  lines.  Evanescent  wave 
components  corresponding  to  the  region  f^+fy  >  1/A^  were  discarded  in  backpropagation 
because  the  propagation  distance  was  much  larger  than  the  wavelength  A.  A  time-delay 
map  was  then  estimated  for  the  backpropagated  wavefront  and  used  to  time-shift  the 
waveforms.  The  subsequent  processing  was  the  same  as  that  for  the  uncompensated  case. 

The  flow-chart  in  Fig.  8  summarizes  the  various  paths  and  steps  taken  to  obtain  images 
from  the  raw  data. 
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III.  METHODS  OF  EVALUATION 

A.  Evaluation  of  Wavefront  Coherence 

Wavefront  coherence  was  measm-ed  by  the  waveform  similarity  factor,  defined  for  N 
waveforms  i  =  1  ~  iV,by 


in  which  the  r^’s  are  adjusted  to  align  all  the  waveforms  in  time  and  to  maximize  the  WSF. 
The  value  of  this  factor  is  easily  shown  to  be  always  between  0  and  1,  with  1  corresponding 
to  the  case  where  all  waveforms  are  exactly  alike  except  for  a  time  delay  and  an  amplitude 
scaling.  The  value  of  the  WSF  depends  on  the  time  delays  Xj  and  will  generally  be  smaller 
if  the  estimated  Xi’s  are  incorrect.  In  the  current  study,  the  Xi’s  were  obtained  using  the 
least-mean-square-error  method,  as  summarized  in  Section  II-A. 


B.  Evaluation  of  Point  Target  Image 

Point  target  images  were  obtained  using  the  single  transmit  image  formation  de¬ 
scribed  in  Section  II. C.  The  point  target  was  illuminated  with  a  focused  transmit  beam 
and  the  corresponding  echoes  received  in  the  two-dimensional  aperture  were  processed  un¬ 
der  different  conditions  of  compensation  to  produce  different  images.  These  images  were 
two-dimensional,  with  one  of  the  dimensions  being  the  array  (azimuthal)  direction  and 
the  other  being  the  time  (range)  direction.  The  evaluation  of  such  images  follows  the 
procedure  we  used  in  a  previous  study  [15].  The  amplitudes  of  the  analytic  envelope  were 
first  projected  using  the  maximum  amplitude  for  each  array  or  time  position,  to  yield  one¬ 
dimensional  amplitude  profiles  along  the  time  or  array  direction.  Effective  widths  were 
then  found  from  these  profiles  at  -10,  -20,  and  —30  dB,  from  the  peak.  The  square 
root  of  the  product  of  the  effective  widths  in  each  direction  was  defined  as  the  effective 
diameter.  The  effective  diameter  was  a  function  of  the  level  from  the  peak,  and  provided 
a  measure  of  the  focal  spot  size  at  different  levels.  Another  measure  was  the  —10  dB  pe¬ 
ripheral  energy  ratio.  This  dimensionless  parameter  was  calculated  by  specifying  an  ellipse 
whose  center  was  at  the  position  of  maximum  amplitude  and  whose  axes  were  half  of  the 
—10  dB  effective  widths  in  each  direction.  The  ellipse  sepeirated  the  image  into  the  central 
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and  the  peripheral  regions,  and  the  ratio  of  the  total  energy  outside  the  ellipse  to  the  total 
energy  inside  the  ellipse  was  the  peripheral  energy  ratio. 

C.  Evaluation  of  Isoplanatic  Patch  Size 

Using  data  measured  as  described  in  Section  IB,  images  of  the  test  point  target . 
were  obtained  at  each  position  after  cross  compensation  using  time-shifts  derived  from  the 
wavefront  measured  with  the  original  point  target.  The  -10  dB  effective  diameter  of  the 
image  was  calculated.  For  self  compensation,  the  compensation  effect  is  the  best,  and  the 
-10  dB  effective  diameter  is  the  smallest.  With  increasing  distance  between  the  test  and 
the  original  point  target,  the  compensation  effect  deteriorates,  and  the  effective  diameter 
increases.  When  the  amount  of  increase  reached  10%,  the  corresponding  distance  between 
the  original  and  the  test  point  targets  was  taken  to  be  half  of  the  isoplanatic  patch  size. 
The  reason  for  the  factor  of  2  is  that  the  isoplanatic  patch  extends  over  both  sides  of  the 
original  point  target. 

Isoplanatic  patch  size  was  evaluated  for  wavefronts  in  the  aperture  as  well  as  after 
backpropagation  by  20  mm.  In  the  latter  case,  all  the  wavefronts  were  backpropagated 
individually  first  ,  and  a  set  of  time  shifts  was  estimated  for  the  backpropagated  wavefront 
corresponding  to  the  original  point  target.  These  time  shifts  were  then  applied  to  the  other 
backpropagated  wavefronts,  which  were  then  focused  and  evaluated  as  before. 

D.  Evaluation  of  Contrast  Resolution 

The  evaluation  of  contrast  resolution  used  b-scan  images  of  the  cyst  phantom  obtained 
with  no  compensation,  with  time-shift  compensation  in  the  aperture,  and  with  backpropa¬ 
gation  followed  by  time-shift  compensation.  The  cyst  phantom  had  a  scatterer-free  cylinder 
that  was  10  mm  in  diameter  and  was  positioned  perpendicular  to  the  transmit  beam.  To 
calculate  the  contrast  ratio,  a  circular  area  of  diameter  6  mm  was  manually  selected  in 
central  part  of  the  cyst  region.  Another  circular  area  of  the  same  diameter  and  12  mm 
above  the  first  circle  was  correspondingly  placed  in  the  region  of  uniform  scattering.  The 
contrast  ratio  was  calculated  as  the  ratio  of  the  average  energy  levels  in  these  two  circles. 
The  area  inside  the  cyst  region  was  manually  selected  because  the  alignment  between  the 
phantom  and  the  scanning  stages  was  not  controlled  exactly.  The  diameter  of  the  circular 
area  was  made  smaller  than  the  known  diameter  of  the  cyst-like  region  in  order  to  exclude 
relatively  large  amplitudes  near  the  edge. 
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IV.  RESULTS 


A  total  of  10  specimens  were  studied  using  the  agar/glass-sphere  phantom.  However, 
three  measurements  were  unusable  due  to  misalignment  of  the  transmit  focus  with  the 
cyst  phantom  or  too  low  SNR’s  caused  by  specimen  attenuation.  Results  obtained  with 
the  remaining  seven  different  specimens  are  presented  here. 

A.  Waveform  Similarity  Factor  and  Time-Delay  Estimation 

Representative  maps  of  estimated  arrival  time  fluctuation  are  shown  in  Fig.  9.  These 
maps  show  essentially  the  same  trend  of  arrival  time  fluctuation  caused  by  the  tissue,  but 
spatial  features  as  well  as  the  numerical  values  of  arrival  time  fluctuations  estimated  using 
the  2"  transducer  data  are  closer  to  those  obtained  with  the  point  target  data  than  those 
in  the  maps  estimated  using  the  1/2"  and  1"  transducer  data. 

The  waveform  similarity  for  the  point  target  and  water  path  propagation  was  0.997, 
very  close  to  the  ideal  value  of  1.0.  The  nonideality  stems  from  the  variation  of  the  response 
characteristics  of  the  elements  in  the  linear  array.  The  waveform  similarity  factor  results 
for  the  seven  specimens  are  tabulated  in  Table  II.  These  results  show  that,  at  the  aperture, 
the  average  waveform  similarity  factor  of  random  scattering  wavefronts  was  0.261  for  1/2" 
transmitter,  0.490  for  1"  transmitter,  and  0.577  for  2"  transmitter.  The  point  target  data 
has  a  much  higher  average  waveform  similarity  factor  at  0.935.  After  backpropagating  the 
wavefronts  by  20  mm,  the  average  waveform  similarity  factor  increased  to  0.275,  0.513,  and 
0.597,  respectively,  for  random  scattering  with  the  1/2",  1",  and  2"  transmitter,  and  to 
0.951  for  the  point  target  data.  These  results  show  that  a  point  target  produces  a  scattered 
wavefront  with  the  highest  spatial  coherence,  but  a  body  of  random  scatterers  illuminated 
with  a  highly  focused  transmit  beam  (such  as  obtained  with  the  2"  transmitter)  can  be 
used  to  produce  a  wavefront  of  sufficient  spatial  coherence  in  the  absence  of  isolated  point 
targets.  These  results  also  show  that  wavefront  backpropagation  removes  some  waveform 
distortion  caused  by  propagation  through  the  inhomogeneous  medium. 

The  time-delay  estimation  results  are  summarized  in  Table  III  in  which  the  rms  value 
of  the  residual  of  fitting  a  spherical  surface  to  the  estimated  arrival  time  surface  is  used 
to  describe  the  fluctuation.  Using  the  results  obtained  with  the  point  target  data  as 
the  standard,  the  rms  arrival  time  fluctuation  estimated  using  random  scattering  was 
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underestimated  by  8.6,  5.7,  and  2.9  ns  for  the  1/2",  1",  and  2"  diameter  transducers, 
respectively.  The  corresponding  average  underestimation  was  reduced  to  6.6,  3.1,  and 
0.4  ns,  respectively,  after  backpropagation  by  20  mm.  These  results,  combined  with  the 
above  waveform  similarity  results,  show  that  a  tight  transmit  focus  is  needed  to  create  a 
highly  coherent  scattered  wavefront  for  the  accurate  measurement  of  time  shifts,  and  that 
backpropagation  increases  the  wavefront  coherence  and  improves  the  accuracy  of  arrival 
time  estimation. 

B.  Point  Target  Images 

Point  target  images  are  shown  in  Fig.  10  to  illustrate  the  compensation  effects  through 
a  tissue  path.  The  uncompensated  image  has  a  large  spread  of  energy  around  the  focal 
point.  The  focus  is  significantly  improved  with  time-shift  compensation  and  the  sidelobe 
energy  is  greatly  reduced.  With  backpropagation  followed  by  time-shift  compensation,  the 
sidelobe  energy  is  further  reduced  but  the  mainlobe  remains  relatively  unchanged.  For 
comparison,  the  corresponding  point  target  image  obtained  through  a  water  path  is  also 
shown.  The  image  obtained  with  backpropagation  followed  by  time-shift  compensation 
approachs  the  image  obtained  through  a  water  path. 

Statistics  are  presented  in  Table  IV  for  the  effective  diameters  of  the  focus  at  differ¬ 
ent  levels  from  the  peak  as  well  as  for  the  -10  dB  peripheral  energy  ratio.  The  results 
show  that  time-shift  compensation  significantly  improves  the  —10  and  —20  dB  effective 
diameter,  but  the  —30  dB  effective  diameter  after  time-shift  compensation  (16.22  mm)  is 
still  much  greater  than  the  corresponding  water-path  value  (7.79  mm).  With  backpropa¬ 
gation  followed  by  time-shift  compensation,  the  effective  diameter  at  -30  dB  is  reduced 
to  8.20  mm,  which  is  very  close  to  the  water-path  value. 

C.  Isoplanatic  Patch  Size 

A  set  of  representative  point  target  images  used  to  determine  the  isoplanatic  patch 
size  is  shown  in  Fig.  11.  These  images  were  obtained  for  different  point  target  positions 
after  time-shift  compensation  using  the  same  time  delay  map  estimated  for  the  first  point 
target  position.  The  deterioration  of  the  compensation  effects  as  the  distance  between 
the  test  point  target  and  the  original  (first)  point  target  increases  is  clearly  visible.  For 
comparison,  the  point  target  image  obtained  without  any  compensation  is  also  shown  in 
the  figure. 
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The  isoplanatic  patch  size  results  are  summarized  in  Table  V.  For  time-shift  compen¬ 
sation  in  the  aperture,  the  average  isoplanatic  patch  size  was  16.7  mm  in  the  azimuthal 
direction  and  39.0  mm  in  the  range  direction.  For  time-shift  compensation  after  backprop- 
agation  by  20  mm,  the  average  isoplanatic  patch  size  increased  to  19.0  mm  in  the  azimuthal 
direction  and  41.4  mm  in  the  range  direction.  This  indicates  that  backpropagation  has 
effectively  extended  the  isoplanatic  patch  size.  In  only  one  case  (npe21)  did  the  isoplanatic 
patch  size  in  the  range  direction  decrease  (from  16.1  mm  to  7.4  mm)  with  backpropagation. 
Further  examination  of  that  case  reveals  that  the  self  compensated  point  target  image  has 
a  —10  dB  focal  diameter  of  2.49  mm  after  backpropagation,  much  smaller  than  the  cor¬ 
responding  value  (3.29  mm)  obtained  with  time-shift  compensation  alone.  Therefore,  the 
threshold  for  determining  the  isoplanatic  patch  size,  which  is  taken  as  1.1  times  the  —10  dB 
focal  diameter  with  self-compensation,  was  correspondingly  lowered.  If  the  same  threshold 
value  were  used  to  determine  the  isoplanatic  patch  size  with  or  without  backpropagation, 
the  isoplanatic  patch  size  would  be  24.2  mm  instead  of  7.4  mm  in  the  range  direction  after 
backpropagation. 

D.  Images  and  Contrast  Ratio 

A  set  of  single  transmit  images  are  shown  in  Fig.  12.  These  images  were  obtained 
with  the  phantom  illuminated  by  the  1/2",  1",  and  2 "diameter  transmitters,  respectively, 
and  the  reception  is  through  a  water  path.  The  data  used  to  form  these  images  were 
from  the  first  scan  position,  so  the  center  of  illumination  is  about  10  mm  away  from  the 
cyst-like  area.  The  reduction  of  beam  size  with  increasing  transmitter  diameter  is  clearly 
demonstrated.  In  the  case  of  the  1/2"  transmitter,  the  illuminated  region  was  so  large 
that  part  of  the  cyst-like  area  appears  in  the  lower  right  corner  of  the  image. 

Another  set  of  single  transmit  images  are  shown  in  Fig.  13.  These  images  are  ob¬ 
tained  with  the  2"  transmitter  through  a  tissue  reception  path  (npe21)  under  different 
conditions  of  beamforming  (UNC,  TSC,  and  BP-I-TSC),  and  through  a  water  reception 
path.  These  images  show  that  the  receiver  beam  was  severely  distorted  by  the  specimen, 
was  significantly  improved  with  time-shift  compensation,  and  was  further  improved  with 
backpropagation  followed  by  time-shift  compensation.  However,  the  results  after  compen¬ 
sation  still  show  more  energy  spread  compared  to  water  path  results. 
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Representative  b-scan  images  that  correspond  to  the  single  transmit  images  shown  in 
Fig.  13  are  shown  in  Fig.  14.  The  contrast  ratio  in  this  example  improved  from  24.0  dB 
to  26.2  dB  using  time-shift  compensation,  and  to  27.5  dB  using  backpropagation  followed 
by  time-shift  compensation.  The  water  path  contrast  ratio  was  much  higher  at  38.0  dB. 
Statistics  of  the  contrast  ratio  results  for  the  seven  tissue  specimens  are  summarized  in 
Table  VI.  The  average  improvement  in  contrast  ratio  is  the  smallest  (2.0  dB)  with  time-shift 
compensation  on  the  1/2"  transducer  data,  and  the  largest  (3.5  dB)  with  backpropagation 
followed  by  time-shift  compensation  on  the  2"  transducer  data. 


22 


VI.  DISCUSSION 


The  reason  for  using  a  configuration  in  which  the  transmit  beam  is  unperturbed  is 
the  need  to  control  the  transmit  beam  size  for  the  purpose  of  this  investigation.  In  in  vivo 
b-scan  imaging,  however,  the  transmit  beam  as  well  as  the  receive  beam  will  be  distorted 
by  the  inhomogeneous  medium,  so  any  compensation  of  the  transmit  beam  must  be  ap¬ 
plied  iteratively  based  on  current  estimation  of  distortion.  To  be  more  specific,  initially, 
a  transmit  focus  could  be  formed  using  geometric  delays,  and  distortion  information  es¬ 
timated  using  the  scattered  wavefront.  Next,  the  distortion  information  can  be  used  to 
improve  the  transmit  focus.  This  would  produce  a  tighter  transmit  focus  for  more  accu¬ 
rate  estimation  of  distortion.  Such  a  process  could  then  be  repeated  until  no  improvement 
is  obtained  or  the  estimated  values  of  distortion  do  not  change  substantially.  To  carry 
out  this  kind  of  iteration,  however,  requires  a  two-dimensional  array  with  electronics  that 
permits  modification  of  time  delay.  If  the  distortion  is  modeled  as  a  time-shifting  screen 
placed  some  distance  away  from  the  aperture,  then  the  transmit  waveforms  should  ideally 
be  individually  adjustable  to  produce  a  predistorted  wavefront  such  that  subsequent  prop¬ 
agation  in  the  inhomogeneous  medium  will  undo  the  distortion  and  a  diffraction-limited 
focus  will  result.  The  predistorted  wavefront  would  be  obtained  by  backpropagating  a 
time-shifted  (but  otherwise  ideally  focused)  wavefront  from  the  position  of  the  hypotheti¬ 
cal  phase  screen  to  the  array.  When  both  transmit  and  receive  beams  are  compensated  for 
distortion,  the  contrast  ratio  improvement  will  be  about  twice  as  much  as  obtained  with 
only  receive  beam  compensation,  because  the  effective  beam  pattern  in  b-scan  imaging 
is  the  product  of  the  transmit  and  receive  beams.  This  implies  that  the  contrast  ratio 
improvements  presented  in  Table  VI  can  potentially  be  doubled  in  practice. 

The  underestimation  of  arrival  time  fluctuation  when  the  transmit  beam  is  broad  (as 
shown  in  Table  III  for  the  1/2"  transducer)  is  noteworthy.  This  phenomenon  has  been 
qualitatively  observed  before  [3].  With  a  broad  beam,  echoes  return  from  a  large  volume 
of  random  scatterers.  Since  the  distortion  effect  for  different  scatterer  positions  is  different 
(with  the  difference  larger  for  scatterers  farther  apart),  the  overall  distortion  effect  for  the 
rf-signal  received  by  a  particular  element  is  in  a  sense  an  average  of  the  different  distortion 
effects.  This  averaging  results  in  the  significant  underestimation  of  arrival  time  fluctuation 
for  a  broad  transmit  beam. 
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A  gradual  variation  of  the  perceived  arrival  time  fluctuation  with  the  point  target 
position  was  observed  in  processing  the  measured  wavefronts  to  derive  the  isoplanatic 
patch  size.  The  variation  consisted  of  changes  in  spatial  feature  details  and  magnitudes 
of  the  fluctuation.  Because  the  relative  position  of  the  receiving  array  and  the  specimen 
remained  fixed  in  our  measurements  for  the  isoplanatic  patch  size,  the  overall  spatial 
features  of  the  arrival  time  fluctuation  remained  unchanged. 

The  size  of  the  isoplanatic  patch  found  in  this  study  implies  that  a  single  distortion 
profile  can  be  used  over  a  reasonable  region  of  about  20  x  40  mm^.  However,  these 
numerical  values  depend  on  many  factors  including  the  severity  of  distortion,  the  distance 
from  the  scattering  region  to  the  aperture,  and  the  threshold  condition  used  to  define  the 
isoplanatic  patch  size.  Examination  of  Tables  V  reveals  that  the  isoplanatic  patch  sizes 
for  datasets  npe21,  npe24,  and  npe27  are  relatively  small.  Correspondingly,  the  waveform 
similarity  factors  for  these  3  sets  of  data  are  low  (Table  II)  and  the  waveform  arrival  time 
fluctuations  are  high  (Table  III).  In  practice,  due  to  the  relatively  slow  speed  of  sound 
propagation  in  tissue,  only  a  limited  number  of  transmissions  can  be  used  for  real-time 
image  formation.  Therefore,  not  many  transmissions  can  be  devoted  to  the  estimation  of 
the  distortion  profile.  To  circumvent  this  difficulty,  a  point  of  interest  may  be  selected, 
and  the  distortion  profile  around  that  point  may  be  estimated  and  used  to  compensate  for 
distortion  in  the  associated  isoplanatic  patch.  This  way  the  frame  rate  is  little  influenced 
and  the  computational  requirements  are  also  reduced. 

The  expansion  of  the  isoplanatic  patch  size  after  backpropagation  merits  further  com¬ 
ment.  If  the  aberration  effect  is  indeed  caused  by  a  time-shifting  screen  placed  at  some 
distance  away  from  the  measurement  aperture,  then,  after  backpropagation  by  the  right 
distance,  the  time  delays  required  to  align  the  waveforms  will  be  the  same  no  matter  where 
the  point  target  is  positioned.  In  this  case,  the  same  set  of  parameters  (backpropagation 
distance  and  arrival  time  fluctuations)  can  be  used  to  correct  for  any  imaging  positions, 
and  the  isoplanatic  patch  size  will  be  very  large.  Therefore,  the  increase  in  the  isoplanatic 
patch  size  after  backpropagation  indicates  that  a  time-shifting  screen  placed  a  certain  dis¬ 
tance  away  is  a  better  model  for  the  aberration  medium.  This  observation  also  indicates 
that  the  isoplanatic  patch  size  is  not  merely  a  function  of  the  aberration  medium  but  is 
also  a  function  of  the  compensation  algorithm. 
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A  comparison,  of  the  images  in  Fig.  12  shows  that  the  beam  pattern  of  the  2"  transducer 
is  only  modestly  narrower  than  that  of  the  1 "  transducer.  A  possible  explanation  is  that 
the  2  "-diameter  transducer  has  a  focusing  lens  that  is  over  1  cm  thicker  at  the  edge  than 
at  the  center.  The  attenuation  through  the  lens  material  presumably  added  a  spatial 
apodization  to  the  transducer  and  effectively  broadened  the  transmit  beam.  The  small 
beamwidth  diflFerence  between  the  1"  and  2"  transducers  is  also  reflected  in  the  very  small 
difference  in  contrast  ratios  (Table  VI)  achieved  with  these  transducers. 

The  results  presented  in  Table  VI  show  that  greater  improvement  is  achieved  with 
a  larger  transmitting  aperture.  An  intuitive  explanation  for  this  is  that  the  scattered 
wavefront  produced  by  the  larger  aperture  is  more  coherent  and  can  be  aligned  with  greater 
accuracy.  Similarly,  wave  backpropagation  removed  some  of  the  amplitude  fluctuation  and 
waveform  distortion  developed  during  propagation  through  the  specimen,  and  this  helped 
to  improve  the  contrast  ratio  further. 

The  compensated  images  still,  however,  have  a  significantly  lower  contrast  ratio  than 
that  of  images  obtained  through  a  water  path.  For  example,  the  average  contrast  ratio 
obtained  with  the  2"  transmitter  after  backpropagation  and  time-shift  compensation  was 
31.3  dB,  whereas  the  water  path  image  in  Fig.  14,  also  obtained  using  the  2"  transmitter, 
had  a  contrast  ratio  of  38  dB.  One  reason  for  the  incomplete  compensation  is  that  the 
backpropagation  model  is  still  an  approximation  of  the  inhomogeneous  medium  and  does 
not  completely  represent  all  the  distortion  effects.  Another  reason  is  the  reduced  signal- 
to-noise  ratio  (SNR)  of  the  received  signals  as  a  result  of  attenuation  in  the  specimen. 
The  average  tissue-path  SNR  obtained  with  the  2"  transducer  was  19  dB  as  measured  by 
digitizing  the  random  signals  from  the  phantom  with  the  pulser  turned  on,  and  by  digitizing 
the  background  noise  with  the  pulser  turned  off.  In  a  water  path,  however,  the  SNR 
obtained  with  the  2"  transducer  (npe284)  was  26.6  dB.  Since  random  noise  is  incoherent 
and  cannot  be  focused,  its  presence  increases  the  baseline  signal  energy  throughout  the 
image,  thereby  reducing  the  image  contrast. 
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VII.  CONCLUSION 


Wavefront  aberration  estimation  and  correction  have  been  investigated  using  random 
scattering  measured  in  a  specially  designed  pulse-echo  configuration  in  which  the  transmit 
beam  was  unperturbed  by  the  aberration  medium. 

From  the  investigation,  the  following  conclusions  can  be  made  about  the  accuracy  of 
aberration  estimation,  the  isoplanatic  patch  size,  and  the  effectiveness  of  correction. 

1)  A  tightly  focused  transmit  beam,  such  as  that  produced  by  the  2"  f/1.5  transducer  in 
this  study,  is  needed  to  produce  a  suflSciently  coherent  scattering  wavefront  so  that  time 
delays  can  be  accurately  estimated. 

2)  Accommodation  of  waveform  shape  distortion  in  time-delay  estimation,  such  as  in  the 
modified  least-mean-square-error  method  that  incorporates  a  global  smoothness  require¬ 
ment  as  described  in  this  paper,  is  needed  to  estimate  the  time-delay  surface  from  a  random 
scattering  wavefront  that  has  propagated  through  distributed  inhomogeneities. 

3)  The  isoplanatic  patch  size  usually  decreases  with  increasing  severity  of  distortion,  as 
shown  in  this  study  for  abdominal  wall  specimens.  The  isoplanatic  patch  size  is  also 
influenced  by  compensation  technique,  as  seen  in  the  increased  isoplanatic  patch  size  for 
backpropagation  followed  by  time-shift  compensation. 

4)  Incorporation  of  waveform  shape  compensation,  such  as  by  using  backpropagation  prior 
to  time-shift  compensation  as  in  this  study,  improves  the  contrast  ratio  over  that  obtained 
by  time-shift  compensation  alone. 

In  the  absence  of  a  priori  knowledge  about  the  propagation  path,  a  tight  transmit 
focus  can  only  be  achieved  by  modifying  iteratively  the  excitation  signals  based  on  previous 
estimation  of  the  aberration.  Ideally,  a  two-dimensional  array  with  programmable  transmit 
waveform  capability  is  needed  for  producing  such  a  tight  focus. 
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APPENDIX 


The  problem  is  to  find  a  spherical  surface  that  best  fits  a  given  arrival  time  surface 
with  an  unknown  offset  tq, 

~  "n)  ~  ^ "k  {Vi  J/o)^  "k  ^0 1 

where  the  sound  speed  c  is  assumed  known,  and  (xo,2/o,-2o)  is  the  unknown  center  of  the 
sphere.  Although  direct  minimization  of  the  sum  of  squared  errors  in  the  above  approxima¬ 
tion  leads  to  a  nonlinear  least-squares  problem,  the  problem  can  be  linearized  by  rewriting 
the  above  equation  as 

c^[T{xi,yi)  -  Tof  ft;  {xi  -  x^f  -t-  {yi  -  yo)^  +  4- 
Rearranging  this  equation  yields 

do  +  XiOi  +  yid2  -f-  nOz  «  bi, 

where 

Ti  =  r{xi,yi), 

9o^  xl  +  y^  +  zl-  {cTof, 
di  =  -2xo, 

^2  =  — 2yo) 
ds  =  2c^ro, 

and 

bi  =  (cTif  -xf  -  y}. 

Therefore,  the  problem  can  be  solved  by  minimizing 

[^0  +  Xidx  -f-  yidz  +  Tidz  —  bi4 

i 

with  respect  to  do,  di,  d2,dz-  This  is  a  linear  problem,  and  xo,  yo,  zq,  and  tq  can  readily  be 
found  from  the  di's. 


Table  I.  Data  Acquisition  Parameters 


Parameter 

.  Transmit  Transducer  diameters 
Transmit  Transducer  focal  length 
2-D  aperture  in  array  direction 
2-D  aperture  in  elevation  direction 
Number  of  scan  lines 
Scan  line  increment 
Signal  center  frequency 
Signal  bandwidth  (—6  dB) 
Sampling  rate 
Sampling  precision 


Value 

0.5",  1",  and  2" 
3" 

66  X  0.72  mm^ 
12  X  1.0  mm2 
20 

1.0  mm 
3.63  MHz 
1.42  MHz 
20  MHz 
12  bits 


Table  II.  Waveform  Similarity  Factor  Statistics.  The  column  under  Pt  are  data  from  a 
point  target,  while  the  columns  under  1/2",  1",  and  2"  are  data  from  random  scattering 
produced  with  transducers  of  the  corresponding  diameters. 


Dataset 


Pt 


In  Aperture 
1/2"  1"  2 


After  Backpropagation 
Pt  1/2"  1"  2 


npel4 

0.977 

0.293 

0.575 

0.715 

0.964 

0.298 

0.579 

0.720 

npel5 

0.956 

0.270 

0.550 

0.674 

0.957 

0.277 

0.560 

0.685 

npel7 

0.979 

0.288 

0.571 

0.733 

0.978 

0.294 

0.575 

0.735 

npe21 

0.863 

0.224 

0.372 

0.389 

0.930 

0.251 

0.442 

0.439 

npe22 

0.935 

0.263 

0.497 

0.544 

0.957 

0.279 

0.535 

0.575 

npe24 

0.893 

0.248 

0.437 

0.474 

0.912 

0.263 

0.449 

0.490 

npe27 

0.941 

0.245 

0.429 

0.511 

0.958 

0.263 

0.449 

0.535 

avg 

0.935 

0.261 

0.490 

0.577 

0.951 

0.275 

0.513 

0.597 

sd 

0.040 

0.023 

0.073 

0.122 

0.021 

0.016 

0.059 

0.109 

Table  III.  Arrival  Time  Fluctuation  Statistics.  The  abbreviations  are  the  same  as  in  Table 
II  while  the  units  are  nanoseconds. 


Dataset 


In  Aperture 
Pt  1/2"  1"  2" 


After  Backpropagation 
Pt  1/2"  1"  T 


npel4 

29.9 

28.6 

28.7 

28.9 

30.6 

27.5 

28.8 

29.3 

npel5 

44.3 

34.9 

40.3 

41.1 

42.4 

35.0 

41.5 

42.0 

npel7 . 

34.8 

30.4 

31.7 

33.1 

35.0 

29.4 

30.8 

32.6 

npe21 

60.9 

48.4 

52.7 

58.0 

61.8 

55.4 

60.6 

65.0 

npe22 

41.6 

30.0 

31.9 

36.0 

41.6 

33.4 

35.6 

39.0 

npe24 

62.9 

52.3 

55.5 

59.4 

61.6 

51.7 

56.5 

60.2 

npe27 

40.2 

30.6 

33.2 

39.5 

42.3 

36.5 

40.0 

44.1 

avg 

45.2 

36.6 

39.5 

42.3 

45.0 

38.4 

42.0 

44.6 

sd 

11.5 

8.9 

9.9 

11.0 

11.3 

10.0 

11.3 

12.4 

31 


Table  IV.  Point-Target  Image  Statistics.  Listed  are  the  average  (avg)  and  standard  devia¬ 
tion  (sd)  of  the  —10,  —20,  and  —30  dB  effective  diameter  in  mm,  and  the  —10  dB  peripheral 
energy  ratio  (PER)  for  7  abdominal  specimens  under  three  conditions:  UNC=uncompensated, 
TSC=time-shift  compensated,  BP-|-TSC=backpropagation  by  20  mm  followed  by  time- 
shift  compensation.  The  corresponding  values  for  a  water  path  are  also  listed. 


UNC  TSC  BP+TSC  Water 


— 10  dB  Dia. 

avg 

sd 

4.38 

3.13 

2.40 

0.08 

2.25 

0.12 

2.25 

—20  dB  Dia. 

avg 

sd 

10.59 

4.43 

4.87 

2.62 

3.70 

0.71 

3.38 

—30  dB  Dia. 

avg 

sd 

18.46 

1.06 

16.22 

2.42 

8.20 

2.68 

7.79 

-10  dB  PER 

avg 

sd 

0.727 

0.537 

0.297 

0.123 

0.203 

0.040 

0.188 

Table  V.  Isoplanatic  Patch  Size  Statistics.  The  units  are  mm. 


p.  In  Aperture  After  Backpropagation 

^  Azimuth  Range  Azimuth  Range 


npel4 

27.4 

56.0 

27.4 

56.0 

npel5 

20.9 

51.5 

23.4 

56.0 

npel7 

27.4 

56.0 

27.4 

56.0 

npe21 

1.6 

16.1 

4.4 

7.4 

npe22 

13.0 

54.1 

15.2 

56.0 

npe24 

16.4 

15.3 

22.2 

19.6 

npe27 

9.9 

24.2 

12.9 

39.1 

avg 

16.7 

39.0 

19.0 

41.4 

sd 

8.7 

18.0 

7.9 

18.9 

Table  VI.  Contrast  Ratio  Statistics.  The  units  are  dB.  Abbreviations:  UNC  =  Uncompen¬ 
sated  data,  TSC  =  time-shift  compensated  data,  BP-fTSC  =  backpropagation  followed  by 
time-shift  compensation,  TSC/2''  delay  =  time-shift  compensation  with  delays  estimated 
using  the  corresponding  2"  transducer  data. 


Dataset 

UNC 

TSC 

1/2" 

BP+ 

TSC 

TSC/ 

2"  delay 

UNC 

1" 

TSC 

BP+ 

TSC 

UNC 

2" 

TSC 

.BP+ 

TSC 

npel4 

24.0 

25.0 

26.1 

24.4 

29.8 

30.6 

30.9 

30.7 

32.9 

33.0 

npel5 

21.9 

22.2 

22.7 

23.7 

28.1 

29.7 

29.6 

29.1 

30.9 

30.5 

npel7 

21.9 

23.8 

24.1 

24.9 

31.0 

31.7 

32.1 

30.5 

33.6 

34.0 

npe21 

13.7 

18.2 

18.0 

18.7 

24.2 

27.6 

28.5 

24.0 

26.2 

27.5 

npe22 

17.9 

19.7 

20.5 

20.6 

28.1 

30.8 

32.4 

28.0 

30.6 

31.9 

npe24 

16.5 

19.7 

19.5 

20.5 

26.8 

30.3 

31.2 

25.0 

29.4 

31.1 

npe27 

20.9 

21.6 

21.5 

20.9 

23.6 

26.1 

26.1 

27.5 

30.3 

31.0 

avg 

19.5 

21.5 

21.8 

22.0 

27.4 

29.5 

30.1 

27.8 

30.6 

31.3 

sd 

3.4 

2.2 

2.6 

2.2 

2.5 

1.8 

2.1 

2.4 

2.2 

1.9 

32 


Fig.  1.  Measurement  Configuration.  A  single-element  transducer  produces  an  unperturbed 
beam  that  illuminates  the  scattering  phantom.  The  phantom  contains  a  scatterer-free 
cylinder  that  emulates  a  blood  vessel  or  a  cyst-like  region.  Echoes  are  received  by  a 
modified  linear  array  after  propagation  through  a  specimen  of  human  abdominal  wall.  The 
transmit  focal  point  is  maintained  in  the  same  spatial  position  by  a  special  transmitter 
mount  (not  shown)  when  the  transmitter  is  changed. 


Fig.  2.  Scan  Movements.  The  linear  array  is  scanned  in  the  elevation  direction  by  scan 
movement  1  to  synthesize  a  2-D  receiving  apertme.  The  transmitter  and  the  linear  array 
(and,  thus,  the  synthesized  2-D  aperture)  are  scanned  together  by  scan  movement  2  to 
form  b-scan  lines. 


Fig.  3.  Diagram  of  Data  Acquisition  System. 
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Fig.  4.  Pulse-Echo  Waveforms  and  Corresponding  Spectra  obtained  through  a  Water  Path 
from  a  Point  Target.  The  top  row  shows  a  waveform  (left)  and  the  corresponding  spectrum 
(right)  from  the  central  array  element  at  the  central  elevation.  The  bottom  row  shows  the 
average  pulse  and  its  spectrum.  The  average  pulse  was  obtained  by  adding  all  the  pulses 
in  the  receiving  aperture  after  the  pulses  have  been  aligned  in  time.  In  the  waveform 
plots,  the  horizontal  axis  is  sample  number  and  the  vertical  axis  is  amplitude  on  a  linear 
scale.  The  transmitter  was  a  2  transducer,  and  the  point  target  was  the  rounded  tip  of  a 
1/8  -diameter  stainless  steel  rod  positioned  so  that  the  tip  echo  had  maximum  amplitude. 
Abbreviations:  pk  =  peak  frequency,  ctr  =  center  frequency  of  —6  dB  passband,  bw  = 
— 6  dB  bandwidth. 
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Fig.  5.  Wavefronts  Measured  through  a  Tissue  Path.  The  first  three  rows  show  wavefronts 
of  scattering  from  the  scattering  phantom  illuminated  with  the  1/2",  1",  and  2"  transducer, 
respectively,  while  the  last  row  shows  the  wavefront  from  the  point  target  illuminated  by 
the  2"  transducer.  The  horizontal  axis  is  time  and  contains  776  samples  obtained  at  a 
20  MHz  sampling  rate.  The  vertical  axis  is  the  array  direction  and  spans  66  elements  with 
a  pitch  of  0.72  mm.  The  rf  signal  amplitudes  are  displayed  using  a  bipolar  logarithmic  gray 
scale  with  40  dB  dynamic  range  for  each  polarity.  The  dashed  lines  indicate  the  segment 
of  signals  used  for  time-delay  estimation  after  geometric  correction. 
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Fig.  6.  Measurements  for  the  Evaluation  of  Isoplanatic  Patch  Size.  The  scattered  wavefront 
from  a  point  target  was  measured  for  an  initial  (original)  position  (indicated  by  solid  lines). 
The  point  target  was  then  moved  in  the  range  or  azimuthal  directions  to  other  (test) 
positions  (indicated  by  dashed  lines)  and  the  scattered  wavefronts  were  measured  again. 


X 


T(AP,0)-T(APoO)  =  At 
PQ=Transmit  Focal  Point 


Fig.  7.  Estimation  of  the  Transmit  Wavefront  Position  along  the  Transmit  Axis  at  Each 
Instant.  The  position  of  the  transmit  focal  point  Po(xo,  yo>  zq)  was  found  by  fitting  a 
spherical  surface  to  the  arrival  time  surface.  For  echoes  arriving  At  later  than  the  echoes 
from  the  focal  point,  the  spatial  origin  is  i{At)  away  from  the  focal  point,  where  i{At) 
satisfies  the  requirement  that  the  propagation  time  along  the  path  APiO  is  At  longer  than 
along  the  path  APqO.  The  angle  a  between  the  transmit  axis  and  the  2:-axis  was  assumed 
known  (35°  for  the  measurements  reported  here). 


Fig.  8.  Flow  of  Data  Processing  and  Image  Formation.  Raw  data  were  first  corrected 
approximately  for  geometric  path  length  difference  using  a  priori  knowledge  about  the 
position  of  transmit  focal  point.  An  arrival  time  surface  was  then  estimated  using  the 
least-mean-square-error  method  with  echoes  arising  from  the  focal  region,  and  a  spherical 
surface  was  fit  to  the  arrival  time  surface  to  determine  the  focal  point  position.  Next,  exact 
geometric  correction  and  dynamic  focusing  were  performed.  A  single  transmit  image  and 
a  b-scan  image  were  then  formed  under  three  different  conditions:  without  compensation 
(UNC),  with  time-shift  compensation  in  the  aperture  (TSC),  and  with  backpropagation 
followed  by  time-shift  compensation  (BP-I-TSC).  The  single  transmit  image  used  aperture 
data  from  a  single  transmit  or  scan  line  position,  while  the  b-scan  image  used  data  from 
different  transmit  or  scan  positions. 


Fig.  9.  Tissue-Path  Arrival  Time  Maps  Estimated  Using  a  Point  Target  and  Using  Random 
Scattering  Produced  by  1/2",  1",  and  2"  Transducers.  The  vertical  coordinate  is  the  array 
direction  and  spans  66x0.72  mm.  The  horizontal  axis  is  the,  elevation  direction  and  spans 
12x1.00  mm.  The  gray  scale  is  linear  with  white  corresponding  to  -flOO  ns  and  black 
corresponding  to  —100  ns.  The  standard  deviation  (indicated  at  the  top  of  each  panel)  and 
the  spatial  features  of  the  arrival  time  map  estimated  from  random  scattering  approach 
those  estimated  using  a  point  target  as  the  transmit  focal  spot  size  gets  smaller.  The 
dataset  was  npe27. 


Fig.  10.  Point  Target  Images  Obtained  through  a  Tissue  Path  (Left  Three  Panels)  and 
through  a  Water  Path  (Rightmost  Panel).  The  tissue  path  images  were  formed  without 
compensation  (UNO),  with  time-shift  compensation  (TSC),  and  with  backpropagation 
followed  by  time-shift  compensation  (BP-I-TSC).  In  each  panel,  the  analytic  envelope  is 
displayed  on  a  40  dB  logarithmic  scale,  and  the  horizontal  (array)  direction  spans  20  mm 
while  the  vertical  (range  or  time)  direction  consists  of  776  time  samples  that  span  29.57  mm. 
The  dataset  was  npe22. 
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Fig.  11.  Point  Target  Images  from  a  Single  Transmit  Position  for  the  Evaluation  of  Iso- 
planatic  Patch  Size.  The  panel  X-Comp[0]  is  the  image  that  results  when  the  wavefront 
measured  at  the  initial  point-target  position  is  compensated  with  time  shifts  derived  from 
the  same  wavefront.  The  panels  X-Comp[2],  [4],  [6],  •  •  [18]  show  images  that  result  when 

wavefronts  measured  with  the  point  target  shifted  in  the  azimuthal  direction  by  1.44  mm, 
2.88  mm,  4.32  mm,  •  •  •,  12.96  mm  are  time-shift  compensated  using  the  time  shifts  derived 
from  X-Comp[0].  The  panel  Uncomp  [0]  shows  the  image  that  results  from  using  the  un¬ 
compensated  wavefront  measured  at  the  initial  point-target  position.  The  scales  are  the 
same  as  in  Fig.  10.  The  dataset  was  npe27. 
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Fig.  12.  Images  for  a  Water  Path  from  a  Single  Transmit  Position  using  Three  Different 
Diameter  Transducers.  For  the  images  from  left  to  right,  the  corresponding  transmitter 
diameter  was  1/2",  1",  and  2".  In  each  panel,  the  analytic  envelope  is  displayed  on  a  30  dB 
logarithmic  scale,  while  the  horizontal  and  vertical  scales  are  the  same  as  in  Fig.  10. 


Fig.  13.  Images  for  a  Water  Path  and  for  a  Tissue-Path  from  a  Single  Transmit  Position 
under  Different  Conditions.  The  2"  transducer  was  used  as  the  transmitter  for  each  image. 
The  left  three  images  were  formed  using  tissue-path  data  without  compensation,  with 
time-shift  compensation  in  the  aperture,  and  with  backpropagation  by  20  mm  followed  by 
time-shift  compensation,  respectively.  The  rightmost  image  was  formed  using  water-path 
data.  The  scales  are  the  same  as  in  Fig.  12.  The  dataset  was  npe21. 


Fig.  14.  B-Scan  Images  Corresponding  to  the  Single  Transmit  Images  in  Fig.  13.  The 
contrast  ratio  (CR)  is  the  difference  between  the  average  energy  levels  for  the  scattering 
region  (el)  and  for  the  scattering-free  region  (e2).  The  scales  are  the  same  as  in  Fig.  10. 
The  dataset  was  npe21. 
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Abstract — Wave  backpropagation  is  a  concept  that  can 
be  used  to  calculate  the  excitation  signals  for  an  array  with 
programmable  transmit  waveforms  to  produce  a  specified 
field  that  has  no  significant  evanescent  wave  components. 
This  concept  can  also  be  used  to  find  the  field  at  a  distance 
away  from  an  aperture  based  on  measurements  made  in 
the  aperture.  For  a  uniform  medium,  three  methods  exist 
for  the  calculation  of  wave  propagation  and  backpropaga¬ 
tion:  the  diffraction  integral  method,  the  angular  spectrum 
method,  and  the  shift-and-add  method.  The  boundary  con¬ 
ditions  that  are  usually  implicitly  assumed  by  these  meth¬ 
ods  are  analyzed,  and  the  relationship  between  these  meth¬ 
ods  are  explored.  The  application  of  the  angular  spectrum 
method  to  other  kinds  of  boundary  conditions  is  discussed, 
as  is  the  relationship  between  wave  backpropagation,  phase 
conjugation,  and  the  time-reversal  mirror.  Wave  backprop¬ 
agation  is  used,  as  an  example,  to  calculate  the  excitation 
signals  for  a  ring  transducer  to  produce  a  specified  pulsatile 
plane  wave  with  a  limited  spatial  extent. 


1.  Introduction 

Three  methods  exist  for  calculating  the  acoustic 
field  in  a  uniform  medium  given  the  vibration  on  a 
planar  surface.  They  are  the  diffraction  integral  method 
[1],  the  angular  spectrum  method  [2], [3],  and  the  more  em¬ 
pirical  shift-and-add  method  [4].  All  these  methods  have 
been  widely  used,  yet  the  underlying  conditions  for  the  ap¬ 
plication  of  each  method  and  the  relationship  between  the 
methods  have  not  been  clearly  addressed.  For  example, 
the  need  for  the  use  of  boundary  conditions  appropriate 
to  the  physical  setting  is  obvious  in  the  diffraction  inte¬ 
gral  method,  yet  the  boundary  conditions  are  not  obvious 
in  the  angular  spectrum  method  and  in  the  shift-and-add 
method.  This  implies  that  some  choice  of  the  boundary 
conditions  has  already  been  made  when  the  latter  two 
methods  are  applied. 

An  application  in  which  the  above  issues  arise  is  the 
determination  of  the  element  excitation  for  an  array  with 
programmable  transmit  waveforms  to  produce  a  specified 
acoustic  wavefront  at  some  position  away  from  the  array. 
This  can  be  accomplished  by  the  backpropagation  of  the 
specified  wavefront  to  the  position  of  each  element  and  the 
assignment  of  the  resulting  signal  as  the  excitation.  The 
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procedure  outlined  here  provides  a  new  approach  to  design 
a  specified  wavefront. 

With  a  planar  array,  backpropagation  can  be  performed 
straightforwardly  using  the  angular  spectrum  method. 
With  a  curved  array  (like  a  convex  or  concave  array),  how¬ 
ever,  the  angular  spectrum  method  is  not  efficient,  and  the 
way  to  backpropagate  with  the  other  two  methods  is  not 
obvious. 

The  purpose  of  this  paper  is  to  show  the  relationship 
between  these  methods  as  well  as  to  answer  the  following 
four  questions. 

1.  What  are  the  implicit  boundary  conditions  in  the 
angular  spectrum  method  and  in  the  shift-and-add 
method? 

2.  How  can  backpropagation  be  accomplished  using  the 
diffraction  integral  method  and  the  shift-and-add 
method? 

3.  Should  the  signals  in  the  shift-and-add  method  be  dif¬ 
ferentiated  with  respect  to  time,  and  should  a  spherical 
spreading  factor  and  an  obliquity  factor  be  included? 

4.  How  may  calculations  of  propagation  and  backpropa¬ 
gation  be  implemented  with  curved  apertures? 

The  computation  of  the  ultrasonic  field  produced  by  a 
transducer  often  uses  the  normal  component  of  particle 
velocity  [5], [6].  This  approach  assumes,  explicitly  or  im¬ 
plicitly,  a  rigid  baffle  boundary  condition,  by  using  a  zero 
value  for  the  normal  velocity  outside  of  the  transducer  sur¬ 
face  in  the  transducer  plane.  Another  common  approach 
is  the  impulse  response  method  [7]  that  makes  the  shift- 
and-add  evaluation  computationally  more  efficient  by  an¬ 
alytically  evaluating  the  diffraction  integral  in  the  time 
domain  for  transducers  of  simple  geometry.  This  approach 
also  assumes  the  rigid  baffle  boundary  condition.  A  recent 
paper  [8]  confirmed  the  equivalence  between  the  impulse 
response  method  and  the  angular  spectrum  method.  How¬ 
ever,  the  use  of  a  rigid  baffle  boundary  condition  may  not 
be  supported  by  experimental  evidence  in  all  cases. 

Measurements  have  been  performed  [9], [10]  to  investi¬ 
gate  the  influence  of  boundary  conditions  on  the  acoustic 
field  produced  by  an  ultrasonic  transducer.  The  usual  con¬ 
figuration  encountered  in  medical  ultrasonic  imaging,  dur¬ 
ing  which  the  transducer  array  is  placed  on  the  subject’s 
skin,  can  be  approximated  by  assuming  the  transducer  to 
be  placed  at  the  surface  of  water.  For  the  latter  situation, 
a  comparison  of  the  water  tank  measurements  with  the¬ 
oretically  computed  results  indicates  that  a  zero  pressure 
boundary  condition  or  a  pressure  release  surface  (p  =  0 
outside  the  transducer  surface)  is  more  appropriate  than 
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a  zero  velocity  boundary  condition  or  a  rigid  baffle  surface 
(u  —  0  outside  the  transducer  surface)  [9]. 

A  rigid  baffle  boundary  condition,  however,  is  more  fre¬ 
quently  used  in  calculating  acoustic  wave  propagation,  in 
spite  of  the  above  noted  experimental  observations.  The 
reasons  for  this  may  include  the  following:  (a)  the  preva¬ 
lent  influence  of  Rayleigh’s  work,  particularly  the  Rayleigh 
integral,  in  acoutics;  (b)  the  applicability  of  a  rigid  baffle 
condition  to  some  common  situations,  such  as  a  speaker 
mounted  on  a  wall;  and  c)  the  small  angles  usually  en¬ 
countered  in  propagation  calculations  that  make  the  use 
of  the  correct  boundary  conditions  relatively  unimportant. 
The  current  paucity  of  data  suggests  the  need  for  further 
detailed  studies,  by  theoretical,  computational,  and  exper¬ 
imental  means,  of  the  fields  produced  by  ultrasonic  arrays 
in  settings  similar  to  those  encountered  in  medical  diag¬ 
nosis.  Such  studies  are,  for  example,  important  for  the  de¬ 
velopment  of  quantitative  ultrasonic  imaging. 

The  remainder  of  this  paper  is  comprised  of  four  sec¬ 
tions:  Theory,  Computations  and  Results,  Discussion,  and 
Conclusion.  In  the  Theory,  the  angular  spectrum  method, 
the  diflFraction  integral  method,  and  the  shift-and-add 
method  are  summarized  in  a  unified  framework  for  cal¬ 
culation  of  wave  propagation  and  backpropagation  in  a 
uniform  medium  in  three-dimensional  or  two-dimensional 
spaces.  An  approximation  is  introduced  to  facilitate  wave 
propagation  calculations  with  a  curved  aperture.  Compu¬ 
tations  based  on  the  analytic  results  are  presented  in  the 
next  section  that  includes,  as  an  example,  a  demonstration 
of  the  use  of  wave  backpropagation  to  obtain  the  excitation 
signals  for  a  ring  transducer  to  produce  a  spatially  lim¬ 
ited  plane  wave.  In  the  Discussion,  the  calculation  of  wave 
propagation  for  other  boundary  conditions  is  considered 
and  the  relationship  between  wave  backpropagation,  phase 
conjugation,  and  the  time-reversal  mirror  is  explored.  The 
main  results  are  summarized  in  the  Conclusion. 


II.  Theory 

Conventions  for  notation  in  this  analysis  are  the  follow¬ 
ing.  Time  dependence  is  assumed  to  be  for  harmonic 
vibrations.  Propagation  or  backpropagation  is  assumed  to 
be  along  the  2;  axis  and  the  variable  2;  is  assumed  to  always 
denote  a  positive  number. 


A.  Angular  Spectrum  Method 

Let  p{x,yyz)  denote  the  complex  temporal  harmonic 
amplitude  of  pressure  at  a  single  frequency  in  a  uniform 
medium.  Given  values  of  p{x^  y,  z)  in  the  xy  plane  at  a 
certain  2:,  the  angular  spectrum  is  defined  [2]  as  the  Fourier 
transform  of  p{x,  y,  z)  with  respect  to  x  and  y.  The  angular 
spectrum  may  thus  be  written 

00  00 

I  I  pix,y,z)e-^^-^-^^+yfyUxdy  .  (1) 

—  00—00 


The  relationship  between  P{fx,  fy>  and  P{fx,  the 

angular  spectrum  at  2  =  0,  is  obtained  from  the  fact  that 
p{x,y,z)  obeys  the  Helmholtz  equation  (V^  -h  k‘^)p  =  0, 
where  k  =  27r/X  is  the  wavenumber.  The  result  is  [2] 

P{fxJy,z)=P{fxJy.O)HiUJy,z)  ,  (2) 

where  is  the  so-called  propagation  function, 

which  can  be  written  explicitly  as 


Hifx,  fy,  z)  = 

/|  +  /2<1/A2 

“  /2  +  f2  >  1/a2  . 


(3) 


In  the  region  +  fy  <  l/A^  the  propagation  function 
has  a  uniform  amplitude  of  1,  and  is  associated  with  what 
are  known  as  homogeneous  waves.  In  the  region  /J  +  /y  > 
1/A^,  the  propagation  function  decays  exponentially  with 
increasing  propagation  distance  2,  and  is  associated  with 
what  are  known  as  evanescent  waves. 

In  arriving  at  (3),  a  choice  was  made  for  the  sign  of  the 
exponent.  The  choice  was  based  on  the  physical  require¬ 
ment  that  homogeneous  waves  with  a  time  dependence  of 
are  advanced  in  phase  to  represent  a  time-delay,  and 
on  the  fact  that  evanescent  waves  decay. 

Equation  (2)  provides  a  way  of  calculating  the  field 
in  an  xy  plane  at  an  arbitrary  2  using  P{fx,fy^O)  or 
p{x,y,0).  Since  the  calculation  in  (2)  is  a  multiplication 
in  the  spatial-frequency  domain,  an  inverse  Fourier  trans¬ 
formation  to  the  space  domain  results  in  a  convolution, 
which  may  be  expressed 


00  00 

p{x,y,z)  =  //  p{xo,  yo,  0)h{x  -  xo,y  -  yo,  z)dxodyo 


where 


00  00 

h{x,y,z)=  I  I  H{^,fy,z)eP<^f^+y^^Uf.dfy 


(4) 


(5) 


is  called  the  impulse  response  function  of  the  propagation 
process.  An  explicit  form  for  h{x,  y,  z)  can  be  obtained  by 
examining  the  angular  spectrum  of  the  free-space  Green’s 
function  g{r)  =  e^^"^  jMxr  where  r  =  \J x^  y^  A-  z^  ^  which 
is  [8], [11] 


G{fx}  fyf  z) 


_ jA _ 

47ry/'l-A2(/|  +  /2) 


(6) 


Since  apparently  dG{fx,  fy,  z)/dz  —  —\H{fx,fy,z),  an 
inverse  Fourier  transform  of  this  relation  with  respect  to 


fx  and  fy  yields 


h{x,y,z) 


2-  ^ 
dz  A inr  ) 


-{'--A-  • 
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This  result  is  the  same  as  that  derived  in  [3]. 
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The  backpropagation  process  is  defined  as  the  calcula¬ 
tion  of  PifxJy.O)  from  P{f:j,,fy,z)\  From  (2),  the  back- 
propagation  function  may  be  defined  as 


g-jfe2^1-A2(/2+/2) 


g-jfc2Vl-A2(/2  +  (/2)^ 
gfez^A2(/2+/2)_l^ 


(8) 


In  the  region  f^  +  fy>  l/A^,  hr(/x,  fy,  -z)  grows  expo- 
nentially  with  2;,  signifying  that  evanescent  waves  get  an 
exponential  amplification  to  compensate  for  the  exponen¬ 
tial  decay  they  undergo  in  forward  propagation.  Because 
H{fx^  fy,  —z)  is  unbounded  in  the  region  >  1/^^? 

however,  its  inverse  Fourier  transform  does  not  exist  so 
that  the  impulse  response  function  of  backpropagation 
cannot  be  defined. 

To  circumvent  this  problem,  consider  two  scenarios  in 
which  backpropagation  might  be  used.  One  is  to  find  ex¬ 
citation  signals  for  the  production  of  a  specified  field.  The 
other  is  to  calculate  the  field  at  some  distance  from  a  mea¬ 
surement  aperture.  At  the  position  of  specification  or  mea¬ 
surement,  because  the  spatial  sampling  rate  is  always  fi¬ 
nite,  the  maximum  spatial  frequency  /max  that  the  data 
can  represent  is  also  finite,  and  the  magnitude  of  angu¬ 
lar  spectrum  is  assumed  to  be  zero  beyond  /max-  In  these 
cases,  therefore,  the  definition  of  H{fx,  fy,—z)  can  be 
modified  by  letting  H{fx,  fy,  -z)  =  0  for  ff  +  >  /^^x- 

Then,  the  inverse  Fourier  transform  H{fx,  fy,  —z)  can  be 
evaluated  at  least  numerically.  If  /max  >  1/A,  using  such 
a  backpropagation  function  still  implies  the  exponential 
amplification  of  the  evanescent  waves  whose  magnitude  of 
spatial  frequency  is  between  1/A  and  /max- 

Evanescent  waves  have  negligible  amplitudes  if  the  dis¬ 
tance  of  propagation  is  more  than  a  few  wavelengths.  For 
example,  for  an  evanescent  wave  component  with  spatial 
frequency  at  2/ A  (double  the  highest  spatial  frequency 
of  nonevanescent  waves),  the  amplitude  after  propagat¬ 
ing  2 A  will  be,  according  to  (3),  reduced  by  = 

3.5  X  10“^®,  which  is  about  —200  dB.  This  huge  attenu¬ 
ation  implies  that  to  produce  evanescent  waves  beyond  a 
few  wavelengths  from  the  aperture  is  almost  impossible, 
and  to  determine  the  existence  of  evanescent  waves  at  a 
few  wavelengths  away  is  also  almost  impossible.  In  this 
case,  the  handling  of  evanescent  waves  becomes  unimpor¬ 
tant,  and  an  arbitrary  function  may  be  used  to  replace 
the  growing  exponential  in  (8)  to  make  the  inverse  Fourier 
transform  well  defined.  One  particular  approximation  that 
leads  to  an  analytic  form  of  the  inverse  Fourier  transform  is 

H{fx,fy,-Z)  =  H*{fx,fy,z) 

^  I ,-ifcVl-A2(/|+/2)^  ^2  +  y2  <  1/^2 

/2  +  /2>1/A2.  ^ 

which  is  obtained  by  replacing  the  growing  exponential 
in  (8)  for  the  evanescent  waves  by  a  decaying  exponential. 
The  corresponding  approximate  backpropagation  impulse 


response  function  is  the  complex  conjugate  of  h(x^y^z), 


h(x,y,-z) 


2  d 

dz  V  47rr 


g-jfcr  /I 

-  -  +  jA:  -  . 

27rr  V r  Jr 


(10) 

The  results  in  (9)  and  (10)  developed  here  for  acoustic  ap¬ 
plications  have  also  been  derived  [12], [13]  independently 
using  similar  reasoning  and  approximation  for  optical  ap¬ 
plications. 

Since  a  forward  propagation  followed  by  a  backpropaga¬ 
tion  of  equal  distance  should  ideally  restore  the  wavefront, 
a  good  way  of  examining  h[x,  ~z)  as  an  approximation 
of  the  ideal  backpropagation  impulse  response  function  is 
to  calculate  the  convolution  of  h{x^y^z)  and  h{x^y^—z) 
with  respect  to  x  and  ?/,  and  to  compare  the  result  with 
the  delta  function,  i.e.,  to  check  the  accuracy  of 


dz  V  47rr  )  dz\  Aixr  ) 


^5{x)5{y)  . 


(11) 


B.  Diffraction  Integral  Method 


Since  p  satisfies  the  Helmholtz  equation  (V^-1-A)^)p  =  0, 
the  value  of  p  at  an  arbitrary  point  f  =  {x,y,z)  can  be 
expressed  as  the  surface  integral  [11] 


pM  =  // 

5o 


^ ^  ^ dp{ro)  ^^^dG{r\fo) 
- p{ro) 


dno 


dno 


dfo 


(12) 

where  Sq  is  an  arbitrary  surface  enclosing  the  point  r,  ro 
is  a  point  on  Sq,  and  no  is  the  outward  pointing  nor¬ 
mal  vector  of  Sq  at  ro-  This  formula  enables  the  calcu¬ 
lation  of  p{f)  from  simultaneous  knowledge  of  both  p{fo) 
and  dp{fo)/dnQ  on  the  surface  Sq.  The  Green’s  function 
G(f|ro)  to  be  used  in  the  above  formula  is  a  general  solu¬ 
tion  of  the  equation  (V^  -h  fc^)G(r|ro)  =  —6{f—  fo),  and 
consists  of  two  terms 


=  9{R)  +  X(10  =  +  X(r)  >  (13) 

where  R  =  |^  —  “rol  is  the  distance  between  f  and  ro, 
g{R)  is  the  free-space  outgoing  Green’s  function  as  already 
defined,  and  x(fO  i®  arbitrary  function  satisfying  the 
Helmholtz  equation  (V^  -h  A^^)x(fO  =  0  in  the  volume  sur¬ 
rounded  by  the  surface  5o.  An  arbitrary  choice  of  x(^  will 
not  change  the  result  computed  by  (12),  but  an  appropri¬ 
ate  choice  of  this  function  eliminates  the  need  of  simul¬ 
taneous  knowledge  of  both  p{fo)  and  dp{ro)/d7io  on  the 
surface  Sq,  thereby  preventing  the  boundary  conditions 
from  being  overspecified. 

Consider  the  case  of  So  being  comprised  of  the  plane 
z  =  0  and  the  hemisphere  at  infinity  enclosing  the  up¬ 
per  half  space.  The  behavior  of  the  pressure  field  at  infin¬ 
ity  is  specified  to  satisfy  the  Sommerfeld  radiation  condi¬ 
tion  [2].  Under  these  conditions,  the  formulae  for  comput¬ 
ing  p{f)  using  either  p{fo)  alone  or  dp{fo)/dno  alone  are 
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well-known  [11], 


or 


z 


,2/0,0) 


—  CO  — oo 


-^dxodyo  , 
(14) 


p{x,  y,  z) 


CO  oo 

^  f  f  dp{xo,yo,zo) 

2ti  J  J  dzo 

—  OO  — OO 


zo=0 


~R~ 


dxodyo  . 


_ (15) 

In  these  expressions,  R  =  \/{x  —  ^o)^  +  iv  ~  Vo)^  +  • 

When  pressure  is  specified  in  the  z  =  0  plane  (known  as 
the  Dirichlet  boundary  condition),  (14)  applies.  A  special 
case  of  the  Dirichlet  boundary  condition  is  the  pressure 
release  surface,  in  which  the  pressure  is  specified  on  the 
transducer  surface  and  is  assumed  to  be  zero  elsewhere  in 
the  transducer  plane.  When  the  normal  derivative  of  pres¬ 
sure  is  specified  in  the  z  =  0  plane  (known  as  the  Neuman 
boundary  condition),  (15)  applies.  A  special  case  of  the 
Neuman  boundary  condition  is  the  rigid  baffle  surface,  in 
which  the  normal  velocity  (which  is  proportional  to  the 
normal  derivative  of  the  pressure)  is  specified  on  the  trans¬ 
ducer  surface  and  is  assumed  to  be  zero  elsewhere  in  the 
transducer  plane.  A  comparison  of  (14)  with  (4)  and  (7) 
immediately  reveals  that  the  angular  spectrum  computa¬ 
tion  corresponds  to  the  Dirichlet  boundary  condition.  This 
conclusion  is  not  surprising,  since  the  angular  spectrum 
method  assumes  that  values  of  p{x,  y,  z)  are  specified  in 
the  plane  z  =  0. 

An  equation  that  describes  backpropagation  and  forms 
a  pair  with  (14)  is  obtained  using  the  approximate  impulse 
response  function  of  backpropagation  in  (10).  The  result  is 


p{xo,yo,  0)  ~  (;^ + J*)  '^^^y  ■ 

—  oo  — oo 

(16) 

The  approximation  in  (16)  is  due,  as  already  noted,  to  the 
exponential  decay  of  evanescent  waves  during  backpropa¬ 
gation  instead  of  exponential  amplification. 


The  inverse  Fourier  transform  of  (17)  with  respect  to  tem¬ 
poral  frequency  then  yields 


p{x,y,z,t)  = 


2itc 


oo  oo 

J  J  p{xo,yo,0,t- R/c)^dxodyo 


(18) 

where  p{x,  y,  z,  t)  denotes  the  temporal  pressure  variation 
at  point  {x,  y,  z)  and  p  denotes  the  differentiation  of  p  with 
respect  to  time.  The  validity  of  the  approximation  leading 
to  (18)  requires  the  relation  X  to  hold  for  all  the  sig¬ 
nificant  frequency  components  contained  in  p{xo,yoj0,t). 
Subsequent  discretization  of  (18)  with  respect  to  {xQ,yo) 
leads  to 


p{x,y,z,t)  =  —  AxoAyo 


^  ^  ^  2/On,  Q,  ^  Rmn/^)  j^2 


(19) 


where  Rmn  =  -  a^om)^  +  [u  -  2/On)^  +  •  Since  (19) 

corresponds  to  time-shifting  and  adding  signals  at  z  =  0 
to  obtain  signals  in  the  plane  at  z,  the  name  shift-and-add 
is  given  to  the  method.  Similarly,  the  spatially  discretized, 
time-domain  version  of  (16)  with  the  assumption  >  A 
can  be  written  as 


z 

p{xo,yo,0,t)  =  Aa;Ay 

x'^^p{xp,yg,z,t  +  Rpg/c)^  >  (20) 

p  q 

where  Rpg  =  y/ {xq  -  Xp)‘^  +  (yo  Vg)'^  T  ^ • 

Both  (19)  and  (20)  come  from  results  that  were  ob¬ 
tained  for  a  planar  geometry.  Thus,  these  equations  also 
only  apply  directly  to  a  planar  geometry,  although  the  re¬ 
quirement  is  limited  to  the  source  aperture  in  the  case 
of  (19)  and  to  the  field  measurement  aperture  in  the  case 
of  (20).  An  approximate  extension  of  these  results  to  more 
general  aperture  shapes  is  possible,  as  described  in  the 
next  section. 


(7.  Propagation  and  Backpropagation  by  Shift-And-Add 

The  shift-and-add  method  can  be  derived  from  the 
diffraction  integrals  by  an  inverse  Fourier  transform  of  (14) 
and  (16)  with  respect  to  time  followed  by  a  discretization 
of  the  integrals  in  space.  An  approximation  is  customarily 
made  under  the  assumption  that  the  distance  R  is  much 
larger  than  the  wavelength  A  (as  is  usually  the  case).  By 
omitting  the  term  1/R  in  the  parentheses  and  denoting 
p{x,  y,  z)  by  pc^{x,  y,  z)  to  signify  that  the  equation  applies 
to  a  single  frequency  component,  (14)  becomes 


D.  Propagation  and  Backpropagation 
in  Two-dimensional  Space 

Assume  now  that  the  field  varies  only  in  the  x  and  z 
directions  and  is  constant  in  the  y  direction  so  the  wave 
propagation  is  in  two  dimensions.  The  angular  spectrum 
method  can  be  applied  to  two-dimensional  wave  propaga¬ 
tion  with  essentially  no  change  in  form.  The  corresponding 
propagation  function  is 


—  jkz  f  f 

Pu>{x,y,z)  ^  J  J  Pu,{xo,yo,  dxodyo 


(17) 


H{f.,z)  = 

^  r  /2  <  i/a2 


(21) 
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> 


and  the  corresponding  approximate  backpropagation  func¬ 
tion  is 

y2< 

J2  >  1/^2  _  ^22) 

Using  an  approach  similar  to  the  one  leading  to  (7),  the 
inverse  Fourier  transform  of  H{fx,  z)  is  found  to  be 


ff(/x 


oo 

h{x,z)=  j  = 

—  OO 


LI. 

2jdz 


rr(l) 

•”o 


where  r  =  y/x^  +  z^  and 

=  MO  +  jNoiO  as  ^  ^  oo 

(24) 

is  the  zeroth-order  Hankel  function  of  the  first  kind.  A 
check  on  the  accuracy  of  the  asymptotic  form  using  the 
IMSL  library  [14]  indicates  that  the  relative  rms  error  for 
the  real  and  the  imaginary  parts  is  less  than  0.1%  for  ^  = 
kr  >  100.  Using  the  asymptotic  form,  the  impulse  response 
function  of  forward  propagation  is  found  to  be 

h(x,z)^^e^^'‘^-^/‘^'>-  .  (25) 

V  Ar 

The  asymptotic  form  of  the  inverse  Fourier  transform  of 
H{fx^  —^)  is  similarly  found  to  be 

h(a;,-z)  .  (26) 

The  above  impulse  response  functions  h{xyz)  and 
h{x,  —z),  as  well  as  those  presented  in  (7)  and  (10)  for  wave 
propagation  in  three-dimensional  space,  were  obtained  as¬ 
suming  a  planar  geometry.  In  practice,  however,  curved 
apertures  may  be  encountered.  To  propagate  from  such  an 
aperture  or  to  backpropagate  fields  measured  by  such  an 
aperture,  an  extension  of  these  results  to  a  curved  geome¬ 
try  is  needed. 

A  rigorous  approach  to  treat  a  curved  aperture  is  to  find 
the  Green’s  function  in  the  form  of  (13)  for  the  particular 
geometry  and  then  to  use  that  Green’s  function  in  (12). 
However,  an  analytic  solution  for  x(f^  is  known  only  for  a 
few  simple  geometries.  Therefore,  this  approach  does  not 
usually  lead  to  uncomplicated  analytic  results. 

Instead,  a  closer  examination  of  the  planar  geometry 
result  using  the  impulse  response  function  proves  worth¬ 
while.  In  the  following,  h{x,  z)  in  (25)  for  wave  propagation 
in  two-dimensional  space  will  be  used  as  an  example  for 
this  examination.  The  field  p(x,  z)  caused  by  a  vibration 
p(x,  0)  at  z  =  0  can  be  expressed  using  h(x,  z), 

oo 

p(x,z)=  J  p{xo,0)h{x  -  xo,z)dxo 

—  OO 
OO 

~  cosOdxo  ,  (27) 

J  vAr 


where  r  =  y^{x  —  Xq)‘^  -h  and  cos0  —  zjv  with  9  being 
the  angle  between  the  normal  of  the  aperture  and  the  line 
that  connects  the  aperture  point  and  the  field  point,  as 
shown  in  Fig.  1.  The  factor  cos  6  is  known  as  the  obliq¬ 
uity  factor.  In  order  to  generalize  this  result  to  an  arbi¬ 
trary  aperture  shape,  the  contribution  of  a  line  segment 
p((ro,  0)^X0  in  the  planar  aperture  to  the  field  at  {x^z)  is 
examined.  Apart  from  terms  related  to  time  delay,  phase 
shifting,  and  amplitude  scaling,  the  term  cos  6  relates  to 
the  directivity  of  the  radiation  pattern  of  the  line  segment. 
Now,  assume  that  the  radiation  pattern  of  a  line  segment 
in  a  curved  aperture  is  little  influenced  by  the  overall  shape 
of  the  aperture  and  is  the  same  as  the  radiation  pattern 
of  a  line  segment  in  a  planar  aperture.  (See  Fig.  1.)  This 
is  probably  a  good  approximation  as  long  as  the  radius  of 
curvature  of  the  aperture  is  much  larger  than  the  wave¬ 
length  (so  that  in  the  vicinity  of  each  line  segment  the 
aperture  is  almost  flat).  Using  this  approximation,  the  con¬ 
tribution  from  individual  line  segments  p(xo,  Z{^)dsQ  can  be 
found  using  the  same  formula  as  in  (27)  and  the  total  field 
is  obtained  by  a  summation.  The  result  is 


oo 

p{x,z)^  j  p{xo,zo)-^e^^’^''~'^^‘^Los9dso  ,  (28) 


where  2;  >  ^o,  r  =  \/{x  —  xq)"^  +  (^  “  dso  is  a  line 
segment  along  the  curved  aperture,  and  9  is  the  angle 
between  the  normal  vector  and  the  vector  r.  A  similar 
approach  was  used  [15]  to  extend  the  impulse  response 
method  (which  assumes  the  rigid  baffle  boundary  condi¬ 
tion)  from  a  planar  aperture  to  a  curved  aperture.  Con¬ 
sider  a  circular  aperture  of  radius  R  that  spans  an  angle  A. 
This  geometry  is  shown  in  Fig.  2.  The  field  produced  by 
this  aperture  can  be  approximately  computed  using  (28), 


A/2 


p{x,  z)  ^  f  p{Rsma^  ~R cos  a) 

J  V  Ar 

-A/2 

,i(kr-n/4)  M+r^-  (X^  +  ZO 


2Rr 


Rda  .  (29) 


Backpropagation  of  a  field  measured  along  a  curve  can 
be  computed  using  the  same  approximation  with  the  obliq¬ 
uity  factor  cos  9  modified  according  to  the  geometry.  The 
formula  is 


00 

p{x,z)^  j  p{xQ,zo)-^e~^^^^~'"/^Los6dsQ  ,  (30) 

—  OO 

where  2:  <  2:0. 

The  time-domain  counterparts  of  (28)  and  (30)  can¬ 
not  be  written  explicitly  because  the  amplitude  scaling 
contains  a  frequency-dependent  factor  1/\/A  =  y/JJc, 
where  /  is  temporal  frequency  and  c  is  sound  speed.  Thus, 
the  computation  of  a  pulse  wavefront  propagating  in  two- 
dimensional  space  needs  to  be  carried  out  for  each  har¬ 
monic  frequency  component  of  the  pulse  individually.  The 
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Plane  Aperture 


cos0dUi 


Fig.  1.  Transition  from  a  planar  aperture  to  a  curved  aperture.  In 
a  planar  aperture  (left),  the  contribution  of  an  aperture  segment 
Uodx  to  a  field  point  in  the  normal  direction  is  denoted  by  dUi.  The 
contribution  to  a  field  point  at  the  same  distance  but  at  an  angle  0 
from  normal  is  given  by  cosOdUi  according  to  (27).  In  a  gradually 
curved  aperture  (right),  the  radiation  directivity  of  cos0  is  assumed 
to  be  little  influenced  by  the  overall  aperture  shape  and  is  used  to 
compute  the  field  contribution  by  a  segment  of  the  aperture. 


z 


Fig.  2.  Geometry  for  the  calculation  of  obliquity  factor  for  a  circu¬ 
lar  aperture.  The  angle  6  is  between  the  normal  direction  and  the 
vector  f  that  goes  from  the  aperture  point  to  the  field  point.  The 
angle  0'  is  between  the  z  axis  and  the  vector  r. 


final  result  in  the  time  domain  can  be  obtained  by  adding 
the  results  for  each  harmonic  frequency  together. 

III.  Computations  and  Results 

Computations  were  made  to  examine  the  validity  of  (11) 
in  a  particular  configuration,  to  investigate  the  effects  of 
the  obliquity  factor  in  three-dimensional  as  well  as  in  two- 
dimensional  wave  propagation,  and  to  demonstrate  the  use 
of  backpropagation  to  determine  element  excitation  of  a 
ring  transducer  for  the  creation  of  a  spatially  limited  plane 
wave. 

A.  Examination  of  the  Validity  of  (11) 

The  computations  were  performed  in  the  following  way. 
First,  complex  functions  h{x^y^z)  and  h{x^y^—z)  were 
evaluated  on  an  (x,  y)  grid  in  the  space  domain.  These  dis¬ 
crete  functions  were  then  Fourier  transformed  using  two- 
dimensional  FFT,  and  the  product  of  the  transforms  was 
then  inverse  Fourier  transformed.  The  use  of  Fourier  trans¬ 
form  in  these  computations  provides  an  efficient  evaluation 
of  the  two-dimensional  convolution  in  (11). 

The  following  parameters  in  the  computation  were  cho¬ 
sen  for  their  applicability  in  medical  ultrasonic  imaging. 
The  center  frequency  was  3.75  MHz.  The  sound  speed  was 
1.5  mm/ (IS.  The  distance  of  propagation  (z)  was  10  mm. 
This  distance  was  intentionally  kept  relatively  short  so 
that  the  impulse  response  function  h{x,y,z)  has  signifi¬ 
cant  values  over  a  relatively  small  area  and  can,  thus,  be 
represented  with  fewer  samples.  This  distance  was  large 
enough,  however,  so  that  evanescent  waves  can  be  safely 
ignored.  An  {x,y)  grid  spacing  of  0.1  mm  x  0.1  mm  was 
used.  This  sampling  interval  is  less  than  half  the  wave¬ 
length,  which  is  0.4  mm  at  3.75  MHz,  and  satisfies  the 
sampling  theorem  requirement. 


Test  computations  were  performed  to  determine  the 
necessary  coverage  of  the  xy  plane,  or  the  grid  size.  The  ex¬ 
pected  angular  spectrum  of  h(x,  y,  z)  has  unit  amplitude 

in  the  region  ^ +  fy  <  1/A  and  decays  exponentially 

beyond  that.  The  graphical  appearance  of  the  magnitude 
of  the  angular  spectrum  is  a  disk.  A  grid  size  of  512  x  512 
points  with  the  grid  spacing  of  0.1  mm  x  0.1  mm  resulted 
in  four  corners  in  the  calculated  angular  spectrum.  These 
were  caused  by  the  truncation  of  h(x,  y,  2;)  in  space  because 
only  an  area  of51.2x51.2  mm^  was  covered  in  this  case 
while  h{x,y,z)  extends  throughout  the  entire  xy  plane. 
When  the  grid  size  was  increased  to  1024  x  1024  points, 
the  four  corners  disappeared,  showing  that  the  spatial  cov¬ 
erage  in  this  case  was  sufficient  to  make  the  spatial  trun¬ 
cation  insignificant.  Zero-padding  to  2048  x  2048  points  in 
space  was  then  used  to  avoid  wraparound  errors  that  may 
occur  in  the  evaluation  of  the  convolution  in  (11)  using 
a  two-dimensional  FFT.  (Wraparound  error  occurs  when 
the  dimensions  of  the  expected  convolution  result  exceed 
the  dimensions  of  the  FFT.) 

The  results  are  illustrated  in  Fig.  3.  These  results  have 
been  normalized  with  respect  to  the  grid  spacing  by  mul¬ 
tiplying  the  FFT  result  and  the  convolution  result  by  spa¬ 
tial  sampling  intervals  Ax  Ay.  The  maximum  magnitude 
of  h{x,  y,  z)  in  panel  (a)  is  at  x  =  y  =  0  and  equals  about 
1/Az,  which  has  the  numerical  value  of  0.25  for  A  =  0.4  mm 
and  z  =  10  mm.  In  panel  (b)  that  shows  the  magnitude 
of  the  angular  spectrum  of  the  impulse  response,  the  di¬ 
ameter  of  the  disk  is  about  1000  FFT  points,  whereas 
2048  FFT  points  corresponds  to  the  spatial  sampling  fre¬ 
quency,  which  is  10  points  per  millimeter.  Thus,  the  ra¬ 
dius  is  at  a  spatial  frequency  of  approximately  2.5  cy/mm, 
which  is  the  expected  1/A.  The  magnitude  in  panel  c)  at 
f^  =  f^  =  0  is  seen  to  be  close  to  1.0,  as  predicted  by  (3). 
Some  ringing  around  the  transition  to  the  evanescent  wave 
region  is  seen  in  the  magnitude.  This  is  a  representation  of 
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Fig.  3.  Computations  associated  with  the  impulse  response  functions  h{x,y,z)  in  (7)  and  h(x,y,  —z)  in  (10).  (a)  Magnitude  of  h{x^y^z)  at 
=  10  mm.  The  function  is  sampled  by  1024  x  1024  points  with  a  spatial  sampling  interval  of  0.1  mm  and  zero-padded  to  2048  X  2048  points, 
(b)  Magnitude  of  the  angular  spectrum  of  y,  z).  (c)  Magnitude  along  the  central  horizontal  line  {fy  =  0)  of  panel  (b).  (d)  Magnitude  of 
the  central  30  X  30  points  that  result  from  convolving  h{x^  y,  z)  with  h{x,  y,  —z),  (e)  The  real  part  of  the  convolution  result  along  the  central 
horizontal  line  (y  =  0) . 


the  Gibbs  phenomenon.  By  applying  an  appropriate,  non- 
rectangular,  window  to  h{x,  y,  z)  in  space  domain,  the  ring¬ 
ing  in  its  angular  spectrum  could  be  minimized  [3],  but  as 
this  was  not  a  concern  in  the  current  investigation,  a  uni¬ 
form  or  rectangular  window  was  used.  The  result  of  convo¬ 
lution  of  h{x,y^  z)  and  h{x^y,  —z)  is  presented  in  panel  d) 
using  the  central  30  x  30  samples  to  show  details  of  the 
result,  since  the  whole  image  of  2048  x  2048  samples  pro¬ 
duces  important  nonzero  values  only  in  the  neighborhood 
of  the  center.  The  integration  of  the  real  part  of  the  convo¬ 
lution  over  the  2048  x  2048  points  was  found  to  be  1.017, 
very  close  to  the  expected  value  of  1.0.  (The  imaginary 
part  is  essentially  zero  other  than  numerical  errors,  due  to 
finite  precision  arithmetic.)  The  distance  between  the  two 
minima  surrounding  the  peak  along  the  central  horizontal 


line,  shown  in  panel  e),  was  found  via  spline  interpolation 
to  be  about  0.66  mm.  This  result  can  be  compared  to  an 
analytical  result,  obtained  by  taking  the  inverse  Fourier 
transform  of  the  disk-shaped  angular  spectrum  (ignoring 
evanescent  waves).  The  result  is  the  well-known  Airy  func¬ 
tion  of  the  form  Ji{kp)/kp  with  p  —  ^/x’^  ^  y^^  and  the 
two  minima  surrounding  the  peak  would  be  1.635 A  apart, 
which  is  0.654  mm  at  the  wavelength  of  0.4  mm.  The  os¬ 
cillatory  tail  in  the  plot  is  due  to  the  sharp  edges  in  the 
spectra  of  h{x,  y,  z)  and  h{x,  y^  —z). 

A  similar  calculation  was  performed  for  —2jkg{x,y^z) 
and  2jkg*{Xyy^z).  These  functions  are  good  approxima¬ 
tions  to  h{x,y^z)  and  h{x^y^—z)  with  the  obliquity  fac¬ 
tor  zjr  removed,  and  would  correspond  to  propagation 
and  backpropagation  by  shift-and-add  using  a  unit  obliq- 
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Magnitudes  along  fy=0 


fx  (cy/mm) 


Spot  Size  0.57  (mm) 


Fig.  4.  Computations  associated  with  the  free-space  Green’s  function. 
The  two  panels  correspond  to  panels  c)  and  e)  in  Fig.  3. 

uity  factor.  The  angular  spectrum  of  -2jkg{x,y,z)  along 
fy  =  0  and  the  result  of  convolution  along  y  =  0,  calcu¬ 
lated  using  the  same  parameters  as  used  for  the  calcula¬ 
tions  in  Fig.  3,  are  presented  in  Fig.  4.  As  can  be  seen 
in  Fig.  4,  the  magnitude  of  the  angular  spectrum,  which 
equals  [1- A2(/|-h/y  increases  from  unity  to  infinity 

as  the  spatial  frequency  increases  from  zero  to  1/A.  The 
singularity  in  the  angular  spectrum  of  g{x^y^z)  is  due  to 
the  uniform  angular  distribution  of  energy  radiated  from 
a  point  source  [16].  The  result  of  convolution  is  more  os¬ 
cillatory  in  Fig.  4  than  in  Fig.  3,  because  the  discontinuity 
in  the  angular  spectrum  at  -h  =  1/A^  is  greater.  The 
integration  of  the  real  part  of  the  convolution  over  the 
2048  X  2048  points  was  found  to  be  1.091,  farther  away 
from  the  expected  1.0  value  than  the  1.017  in  Fig.  3.  The 
distance  between  the  two  minima  surrounding  the  peak  is 
found  to  be  about  0.57  mm,  which  is  14%  smaller  than  the 
spot  size  of  Fig.  3. 

B.  Use  of  Backpropagation  to  Determine  Element 
Excitation  for  the  Creation  of  a  Specified  Field 

The  problem  that  initially  stimulated  the  investigation 
presented  here  is  to  prescribe  a  set  of  waveforms  that  may 
be  applied  to  a  segment  of  the  elements  of  a  ring  transducer 
to  create  a  spatially  limited  plane  wave  at  the  center  of  the 
ring.  This  ring  transducer  was  developed  by  a  team  at  Uni¬ 
versity  of  Rochester  for  ultrasonic  scattering  and  imaging 
studies.  The  ring  transducer  has  a  diameter  of  150  mm 
and  consists  of  2048  elements,  each  able  to  transmit  and 
receive  independently  with  a  center  frequency  of  2.43  MHz 
and  a  -6  dB  bandwidth  of  1.71  MHz.  The  transmit  wave¬ 
form  (pulse  shape)  is  independently  specifiable  for  each 
element,  a  unique  capability  not  found  in  commercial  ul¬ 
trasonic  imaging  systems.  A  diagram  of  the  geometry  is 
shown  in  Fig.  5. 


0.23  mm 


135®  (768  elements) 


Fig.  5.  Ring  transducer  and  spatially  limited  plane  wave.  The  ring 
is  composed  of  2048  elements,  with  the  longer  dimension  of  each  el¬ 
ement  being  parallel  to  the  ring  axis  or  perpendicular  to  the  paper 
plane.  A  wavefront  whose  amplitude  a{x)  along  the  diameter  is  uni¬ 
form  at  the  central  region  of  50  mm  is  to  be  produced.  A  cosine 
roll-off  extending  1  mm  on  each  side  of  the  central  region  is  used  to 
provide  a  smooth  transition  in  amplitude. 

The  field  along  the  150  mm  diameter  was  discretized 
into  1000  points.  This  discretization  meets  the  sampling 
requirement  for  waves  with  a  frequency  up  to  5  MHz.  The 
field  was  specified  to  have  a  uniform  amplitude  for  the 
central  50  mm  with  an  additional  1  mm  cosine  roll-off  at 
each  end,  and  to  have  zero  amplitude  along  the  remain¬ 
der  of  the  diameter.  The  pulse  waveform  was  specified  to 
have  a  Gaussian  envelope,  a  center  frequency  of  2.43  MHz, 
and  a  —6  dB  bandwidth  of  1.71  MHz.  The  field  was  thus 
represented  as 

p{x,  t)  =  a{x)e~^^ sin(27r/ot)  ,  (31) 

where  a{x)  is  as  shown  in  Fig.  5,  fo  =  2.43  MHz,  and  ar 
was  adjusted  to  meet  the  -6  dB  bandwidth  requirement. 
The  pulse  wavefront  was  temporally  sampled  at  20  MHz. 

The  basic  approach  was  to  backpropagate  the  specified 
field  to  the  position  of  each  element  and  prescribe  the  cal¬ 
culated  vibration  as  the  excitation  signal  for  that  element. 
The  field  produced  by  the  prescribed  excitation  was  then 
calculated  (using  forward  propagation)  and  compared  to 
the  specified  field.  Experimental  verification  is  beyond  the 
scope  of  this  paper  (but  is  planned  for  inclusion  in  another 
manuscript) . 

In  the  computations,  two-dimensional  wave  propaga¬ 
tion  was  used  as  opposed  to  three-dimensional  wave  prop¬ 
agation  for  three  reasons.  First,  putting  the  physical  ring 
transducer  aside,  a  check  of  the  approximations  used  to  ob¬ 
tain  (28)  and  (30),  i.e.,  to  see  how  well  a  backpropagation 
to  a  curved  aperture  followed  by  a  forward  propagation 
from  that  curved  aperture  restores  the  original  wavefront, 
is  of  interest.  Second,  the  far-field  point  of  each  element  in 
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the  elevation  (height)  direction,  estimated  using  in 

which  h  =  25  mm  is  element  height  and  A  =  0.6  mm  is  the 
center  frequency  wavelength,  is  260  mm  away  from  the 
element.  This  is  significantly  larger  than  the  ring  radius 
(75  mm)  and  indicates  that  each  element  is  closer  to  a  line 
source  than  a  point  source  and  that  wave  spread  in  the  ele¬ 
vation  direction  is  small  within  the  range  of  consideration. 
Third,  the  ring  transducer  is  essentially  a  one- dimensional 
array  so  that  each  element  of  the  transducer  can  only  act 
as  a  line  transmitter  or  a  line  receiver  (of  finite  length) 
and  cannot  be  divided  into  smaller  units  along  the  eleva¬ 
tion  direction.  Thus,  three-dimensional  wave  propagation 
results  cannot  be  applied  directly  to  the  ring  transducer 
system. 

The  results  of  the  two-dimensional  wave  propagation 
and  backpropagation  calculations  described  above  are 
shown  in  Fig.  6.  The  number  of  elements  used  for  the  pro¬ 
duction  of  the  spatially-limited  plane  wave  was  768  (out 
of  the  total  of  2048).  This  number  was  chosen  based  on 
tests  that  showed  that  the  calculated  excitation  falls  off 
to  a  negligible  amplitude  beyond  that  aperture  size.  To 
reduce  computation,  the  forward  propagation  of  the  exci¬ 
tation  signals  was  made  with  a  coarser  spatial  sampling, 
in  which  100  positions  were  employed  to  cover  the  ring 
diameter  (150  mm).  A  small  error  (about  0.9%  in  magni¬ 
tude)  was  observed  at  the  border  of  the  spatially  limited 
plane  wavefront.  The  error  is  attributed  to  the  existence 
of  a  nonnegligible  amount  of  evanescent  wave  energy  in 
the  specified  field.  In  fact,  when  the  width  of  cosine  roll¬ 
offs  was  increased  from  1  mm  to  5  mm,  the  error  reduced 
dramatically  to  0.03%. 

In  the  above  computation,  (29)  was  used  to  calculate 
the  forward  propagation  from  the  circular  aperture.  For 
comparison,  computations  were  also  made  to  investigate 
the  error  that  can  be  caused  by  using  incorrect  obliquity 
factors.  Two  incorrect  obliquity  factors  were  considered: 
unity  and  cos  Q’  with  &  being  the  angle  between  the  2:  axis 
and  the  vector  r  (Fig.  2).  In  Fig.  7,  the  amplitudes  across 
the  diameter  are  plotted  for  the  originally  specified  field 
and  for  three  other  fields,  all  obtained  with  backpropaga¬ 
tion  of  the  specified  field  to  the  ring  followed  by  forward 
propagation,  but  using  the  three  different  obliquity  factors 
during  forward  propagation.  The  error  caused  by  the  use 
of  incorrect  obliquity  factors  is  clearly  shown. 


IV.  Discussion 

A.  Nonplanar  Geometry  and  Efficiency  of  Computation 

The  extension  of  wave  propagation  formulae  from  a  pla¬ 
nar  aperture  to  curved  apertures  was  only  presented  for 
two-dimensional  space.  However,  the  same  kind  of  exten¬ 
sion  applies  as  well  to  calculations  of  three-dimensional 
wave  propagation,  using  either  the  diffraction  integral 
method  or  the  shift-and-add  method,  but  not  the  angular 
spectrum  method  (because  the  last  method  is  inherently 
based  on  a  planar  geometry).  The  extension  is,  however, 


approximate  and  applies  only  to  curved  apertures  whose 
radius  of  curvature  is  much  greater  than  the  wavelength. 

A  different  approach  is  possible  for  wave  propagation 
calculations  associated  with  a  nonplanar  geometry.  The 
approach  initially  propagates  the  field  from  a  curved  sur¬ 
face  to  a  close-by  planar  surface.  As  long  as  the  distance  of 
propagation  is  small,  a  simple  time  delay  or  a  phase  shift 
(for  a  single  frequency  component)  can  be  used  to  perform 
the  propagation  approximately.  Then,  the  propagation  or 
backpropagation  starting  from  the  planar  surface  can  be 
performed  exactly  as  described  here.  This  approach  is  sim¬ 
ilar  to  the  thin  lens  approximation  used  in  optics. 

The  numerical  efficiency  of  the  angular  spectrum 
method  for  planar  geometry  deserves  emphasis.  For  each 
single  frequency,  the  propagation  from  one  plane  to  an¬ 
other  plane  requires  only  a  Fourier  transform,  a  com¬ 
plex  multiplication,  and  an  inverse  Fourier  transform. 
The  Fourier  transforms  can  be  implemented  using  a  Fast 
Fourier  Transform  algorithm.  As  a  result,  for  the  cal¬ 
culation  of  propagation  or  backpropagation  between  two 
planes,  the  angular  spectrum  method  is  computationally 
the  most  efficient.  However,  this  method  is  not  applicable  if 
the  propagation  or  backpropagation  starts  from  a  curved 
boundary,  in  which  case  the  other  two  methods  may  be 
employed  with  the  above-mentioned  extension. 


B.  Angular  Spectrum  Method  with 
Other  Boundary  Conditions 


The  angular  spectrum  calculation  in  the  theory  sec¬ 
tion  for  three-dimensional  wave  propagation  proceeds  un¬ 
der  the  assumption  that  values  of  pressure  p{x,y,z)  are 
given  in  the  aperture.  The  calculation  requires  modifica¬ 
tion  when  other  kinds  of  boundary  conditions  are  specified. 
Consider  a  planar  surface  of  local  reaction  [11]  with  a  uni¬ 
form  impedance  Z,  The  boundary  condition  imposed  by 
the  surface  is 

tp(x,  y,  0)  +  (x,  y,  0)  =  b{x,  y)  ,  (32) 

where  Vz  is  the  normal  velocity  which  is  related  to  the 
normal  derivative  of  pressure  through 


Vz{x,y,0) 


1  dp{x,y,z) 
jiop  dz 


(33) 


and  b{xy  y)  is  the  boundary  value  that  is  usually  nonzero 
on  the  surface  of  the  transducer  and  zero  elsewhere.  Sub¬ 
stituting  (33)  into  (32),  taking  Fourier  transform  with  re¬ 
spect  to  X  and  y^  and  using  the  fact  that  the  Fourier 
transform  of  dpjdz  at  2:  =  0  equals  P{fx^fy^^)  times 

jk^l  -  A2(/2  +  /2),  as  is  evident  from  (2),  yields 


P(/../„0)  =  T - =  ,  (34) 

where  B{fxyfy)  is  the  Fourier  transform  of  5(x,y).  Once 
P{fx,  fyi  0)  is  available,  P{fx:  fy>'P)  s^t  an  arbitrary  2:  can 
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Fig.  6.  Wavefront  backpropagation  applied  to  calculate  ring  transducer  excitations  for  the  production  of  a  spatially  limited  plane  wave. 
Prom  left  to  right  and  top  to  bottom,  respectively,  are;  the  specified  field,  the  calculated  excitation  signals,  the  field  obtained  by  forward 
propagating  the  excitation  wavefront,  and  the  difference  between  the  calculated  and  the  specified  field.  In  all  the  panels,  horizontal 
time  sampled  at  20  MHz.  In  the  top-right  panel,  the  vertical  axis  corresponds  to  element  number  and  the  display  is  in  log  scale  with  a  50  dB 
range  for  each  polarity.  In  the  other  panels,  the  vertical  axis  corresponds  to  100  spatial  positions  equally  spaced  along  the  ring  diameter 
and  the  display  is  on  a  linear  scale  with  gray  being  zero,  white  and  black  being  positive  and  negative  maximum  amplitudes,  respectively. 
The  amplitude  range  is  indicated  in  the  lower  right  corner  of  each  panel. 


be  obtained  using  (2).  When  Z  ^  0,  the  above  result 
approaches  that  for  a  pressure  release  surface,  while  when 
Z  —>  00,  the  above  result  approaches  that  for  a  rigid  baffle 
surface. 

In  the  particular  case  of  Z  oo,  B{fx,fy)  becomes 
VzifxJy),  and  (34)  becomes 


Pif.Jy.O) 


pO'^z  {fxt  fy) 

^l-A2(/2  +  /2) 


(35) 


This  equation  shows  that  a  normal  velocity  boundary  con¬ 
dition  Vz{x,  y,  0)  can  be  effectively  transformed  into  a  pres¬ 
sure  boundary  condition  p{x,y,0).  Multiplying  (35)  with 
H{fx,fy,z)  in  (3),  performing  inverse  Fourier  transform 


with  respect  to  fx  and  fy,  and  using  the  angular  spectrum 
of  free-space  Green’s  function  in  (6),  yields 


p{x,y,z)  =  j  j 


Vz{xo,yo,0)e^’^^ 


R 


which  is  the  well-known  Rayleigh  integral 
domain. 


■dxodyo  ,  (36) 

in  the  frequency 


C.  Backpropagation,  Phase  Conjugation, 
and  the  Time-Reversal  Mirror 

The  relationship  between  wave  backpropagation,  phase 
conjugation,  and  the  so-called  time-reversal  mirror  merits 
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Fig.  7.  Comparison  of  using  correct  and  incorrect  obliquity  factors 
during  forward  propagation  of  excitation  wavefront  from  the  ring. 
Prom  left  to  right  and  top  to  bottom,  respectively,  are:  the  specified 
field  amplitude  and  the  calculated  field  amplitudes  using  the  correct 
obliquity  factor  cos^  as  defined  in  Fig.  2,  using  a  unity  obliquity 
factor,  and  using  the  obliquity  factor  cos  9'  also  as  defined  in  Fig.  2. 
The  horizontal  axis  in  these  plots  corresponds  to  100  positions  along 
the  ring  diameter. 


comment,  since  these  techniques  are  becoming  more  widely 
used  and  the  analysis  of  wave  propagation  and  backprop- 
agation  in  this  paper  provides  an  opportunity  to  eluci¬ 
date  the  similarity  and  dissimilarity  of  these  techniques. 
For  the  purpose  of  discussion,  assume  a  real  pressure  field 
p(x,  y,  2:,  t)  is  given  in  an  xy  plane,  and  its  temporal  Fourier 
transform  is  denoted  by  Puj{x^  y,  z).  Because  p{x,  y,  z,  t)  is 
a  real  quantity,  the  result  of  phase  conjugation  p*  (x,  p,  z) 
equals  the  Fourier  transform  of  p{x,y,z,—t).  Therefore, 
aside  from  differences  in  implementation  [17],  phase  con¬ 
jugation  and  time-reversal  differ  only  in  that  the  former 
operates  on  a  single  frequency  component  whereas  the  lat¬ 
ter  is  the  time-domain  representation  of  the  same  opera¬ 
tion  on  all  the  frequency  components.  However,  the  time- 
domain  approach  has  extra  flexibility:  the  signal  can  be 
temporally  windowed  so  that  only  a  selected  part  of  the 
signal  is  time-reversed.  The  application  of  phase  conju¬ 
gation  or  the  time-reversal  mirror  usually  consists  of  two 
stages.  First  a  field  is  measured.  Then  the  measurement  is 
either  phase  conjugated  or  time-reversed  and  retransmit¬ 
ted.  The  retransmission  involves  a  physical  propagation  in 
the  medium  so  the  medium  must  be  present.  In  contrast, 
wave  backpropagation  is  a  computational  process  in  which 
no  physical  propagation  takes  place,  so  the  medium  need 
not  be  present. 

To  examine  further  the  relationship  between  backprop¬ 
agation  and  phase  conjugation  (or  the  time-reversal  mir¬ 
ror),  consider  a  harmonic  wavefront  with  complex  ampli¬ 
tude  p{x,  y,  z)  at  a  given  z.  The  result  of  approximate  back- 
propagation  (in  which  evanescent  waves  decay  instead  of 
amplify  with  backpropagation)  of  p(x,  y,  z)  to  the  plane  z  = 


0  can  be  computed  using  H*  given  in  (9).  The  result  is 

pBp{f.Jy.O)  =  P{f,Jy,z)H%f^Jy,z)  ,  (37) 

where  the  subscript  BP  denotes  backpropagation.  On  the 
other  hand,  the  same  wavefront  may  be  phase  conjugated 
and  then  retransmitted  to  produce  a  wave  that  propagates 
back  to  the  plane  at  2;  =  0.  In  this  case,  although  the  phys¬ 
ical  propagation  is  in  the  —2:  direction,  the  propagated  ho¬ 
mogeneous  waves  {fxPfy  ^  with  a  time  dependence 
of  are  advanced  in  phase  to  represent  a  time  delay, 

and  the  amplitude  of  the  propagated  evanescent  waves  de¬ 
cay  [18].  Therefore,  the  propagation  function  is  exactly  the 
same  as  H{fx^  fy^  z)  given  by  (3).  Since  the  angular  spec¬ 
trum  of  the  phase-conjugated  wavefront  p'^{x^y,z)  is 

Ppc{f.,fy,z)^  1 1  p*ix,y,z)e-^^<-f^+yfykxdy 

=  P*i-f.,-fy,z)  ,  (38) 

where  the  subscript  PC  denotes  phase  conjugation,  the 
result  of  propagating  a  phase-conjugated  wavefront  in  the 
—2:  direction  is 


Ppcifxjy,0)  =  P*{-U,-fy,z)Hif,,fy,z) 

^[Pi-f.,-fy.z)H*i-f,,-fy,z)r 

—  ^Bp(“/x)  0)  •  (39) 

Interpreted  in  the  space  domain,  this  result  shows  that 

Ppc (x,  y,  0)  =  p^p (x,  y,  0)  ,  (40) 

which,  in  turn,  shows  that  the  time-domain  result  of  prop¬ 
agating  a  phase-conjugated  wavefront  is  the  time-reversal 
of  the  approximately  backpropagated  wavefront. 

A  wavefront  p{x,  y,  0)  that  starts  at  plane  2:  =  0,  prop¬ 
agates  to  the  plane  at  2:  and  is  measured  there  can,  there¬ 
fore,  be  fully  recovered  by  backpropagation  to  2:  =  0  if  ex¬ 
ponential  amplification  of  evanescent  waves  is  performed 
during  backpropagation.  The  backpropagated  wavefront 
can  have  a  resolution  higher  than  that  limited  by  H- 
fy  <  1/A^,  as  demonstrated  in  nearfield  acoustic  holog¬ 
raphy  where  measurements  were  used  to  reconstruct  the 
field  [19].  On  the  other  hand,  the  p(a;,y,  2;)  wavefront  can 
be  phase  conjugated  and  retransmitted  to  propagate  back 
to  the  plane  2:  =  0.  In  the  latter  case,  the  wavefront  re¬ 
sulting  at  2:  =  0  is  approximately  the  time-reversal  of  the 
initial  wavefront.  The  approximation  results  because  the 
evanescent  waves  have  decayed  twice  during  this  process: 
once  from  0  to  2:  and  once  from  2:  to  0.  Because  of  the  de¬ 
cay,  the  field  obtained  by  retransmission  has  a  wavelength- 
limited  resolution. 

A  matched  filter  model  was  employed  in  [17]  to  explain 
that  a  time-reversal  mirror  can  be  used  to  achieve  a  point 
focus  by  stating  that  the  contribution  from  each  element 
is  maximized  at  the  same  instant  of  time  at  a  certain  field 
point.  However,  this  statement  applies  to  any  field  point 
within  a  distribution  of  random  scatterers  and  does  not  ex¬ 
plicitly  identify  the  need  for  an  isolated  point  target.  The 
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need  for  a  point  target  also  limits  the  focus  produced  by 
time-reversal  to  be  at  the  position  of  the  point  target,  and 
no  way  apparently  exists  to  steer  the  transmit  beam  ef¬ 
fectively  to  other  locations  (which  is  essential  for  imaging 
purpose)  unless  the  inhomogeneous  medium  is  modelled 
such  as  by  a  phase  screen  placed  some  distance  away  from 
the  aperture  [20].  Another  impression  obtained  from  the 
matched-filter  explanation  is  that  the  result  of  retrans¬ 
mission  contains  a  temporal  convolution  of  the  form  of 
hi{ro,T  -  t)  *  hi{fo,t),  where  h^{fo,t)  is  the  acoustic  im¬ 
pulse  response  from  field  point  fo  to  transducer  element  i. 
See  (9)  of  [17].  Therefore,  pulses  may  become  longer  even 
if  the  acousto-electric  impulse  response  of  the  elements 
hf^{t)  is  an  ideal  6{t).  The  analysis  presented  above  shows 
that  this  does  not  happen  in  a  uniform  medium  if  evanes¬ 
cent  waves  are  ignored. 

To  summarize,  backpropagation  differs  conceptually 
from  phase  conjugation/time- reversal  mirror.  Backprop¬ 
agation  is  a  computational  process  and  can  reconstruct  a 
field  with  a  resolution  higher  than  that  limited  by  wave¬ 
length,  whereas  phase  conjugation/time-reversal  mirror  is 
a  physical  process  and  the  retransmitted  field  always  has 
a  resolution  that  is  wavelength  limited.  When  applied  for 
imaging  through  an  inhomogeneous  medium,  backprop¬ 
agation  assumes  knowledge  or  requires  a  model  of  the 
medium,  and  can  be  applied  for  beam  steering  and  imag¬ 
ing,  whereas  phase  conjugation  or  time  reversal  requires 
that  the  medium  be  present  and  that  an  isolated  point 
target  exist  to  form  a  focus  only  at  that  point. 


V.  Conclusion 

The  computation  of  wave  propagation  and  backprop¬ 
agation  in  a  three-dimensional  or  a  two-dimensional  uni¬ 
form  medium  has  been  considered.  Three  methods,  ie., 
the  angular  spectrum  method,  the  diffraction  integral 
method,  and  the  shift-and-add  method,  have  been  exam¬ 
ined  and  their  relationships  explored.  The  straightforward 
application  of  the  angular  spectrum  method  was  shown 
to  correspond  to  Dirichlet  boundary  condition.  The  shift- 
and-add  method  was  shown  to  require  a  temporal  differ¬ 
entiation  of  the  signal  as  well  as  amplitude  scaling  and 
multiplication  with  an  obliquity  factor  in  order  for  the  re¬ 
sults  to  be  the  same  as  obtained  with  the  other  two  meth¬ 
ods.  By  replacing  the  exponential  amplification  of  evanes¬ 
cent  waves  during  backpropagation  with  an  exponential 
decay,  results  were  obtained  for  backpropagation  using  the 
diffraction  integral  method  and  the  shift-and-add  method 
in  addition  to  the  angular  spectrum  method.  Wave  prop¬ 
agation  in  two-dimensional  space  was  shown  to  differ  in 
several  aspects  from  three-dimensional  wave  propagation. 
In  particular,  no  corresponding  formula  for  shift-and-add 
method  exists  in  the  two-dimensional  case.  The  diffraction 
integral  method  for  two-dimensional  propagation  has  been 
extended  to  apertures  of  arbitrary  shape  by  changing  the 
obliquity  factor  according  to  the  geometry.  This  extension 
is  based  on  the  assumption  that  the  aperture  is  smooth 


and  that  the  contribution  of  each  element  on  the  aperture 
surface  to  a  field  point  is  little  influenced  by  the  overall 
shape  of  the  aperture.  Computations  have  been  used  to 
illustrate  the  above  results. 
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An  inverse  scattering  method  that  uses  eigenfunctions  of  the  scattering  operator  is  presented.  This 
approach  provides  a  unified  framework  that  encompasses  eigenfunction  methods  of  focusing  and 
quantitative  image  reconstruction  in  arbitrary  media.  Scattered  acoustic  fields  are  described  using  a 
compact,  normal  operator.  The  eigenfunctions  of  this  operator  are  shown  to  correspond  to  the 
far-field  patterns  of  source  distributions  that  are  directly  proportional  to  the  position-dependent 
contrast  of  a  scattering  object.  Conversely,  the  eigenfunctions  of  the  scattering  operator  specify 
incident- wave  patterns  that  focus  on  these  effective  source  distributions.  These  focusing  properties 
are  employed  in  a  new  inverse  scattering  method  that  represents  unknown  scattering  media  using 
products  of  numerically  calculated  fields  of  eigenfunctions.  A  regularized  solution  to  the  nonlinear 
inverse  scattering  problem  is  shown  to  result  from  combinations  of  these  products,  so  that  the 
products  comprise  a  natural  basis  for  efficient  and  accurate  reconstructions  of  unknown 
inhomogeneities.  The  corresponding  linearized  problem  is  solved  analytically,  resulting  in  a  simple 
formula  for  the  low-pass-filtered  scattering  potential.  The  linear  formula  is  analytically  equivalent  to 
known  filtered-backpropagation  formulas  for  Born  inversion,  and,  at  least  in  the  case  of  small 
scattering  objects,  has  advantages  of  computational  simplicity  and  efficiency.  A  similarly  efficient 
and  simple  formula  is  derived  for  the  nonlinear  problem  in  which  the  total  acoustic  pressure  can  be 
determined  based  on  an  estimate  of  the  medium.  Computational  results  illustrate  focusing  of 
eigenfunctions  on  discrete  and  distributed  scattering  media,  quantitative  imaging  of  inhomogeneous 
media  using  products  of  retransmitted  eigenfunctions,  inverse  scattering  in  an  inhomogeneous 
background  medium,  and  reconstructions  for  data  corrupted  by  noise.  ©  1997  Acoustical  Society 
of  America.  [80001-4966(97)02308-4] 
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INTRODUCTION 

This  paper  presents  a  new  inverse  scattering  method  that 
employs  the  focusing  properties  of  certain  acoustic  fields  ob¬ 
tained  by  retransmitting  eigenfunctions  of  the  scattering  op¬ 
erator. 

Eigensystem  decomposition  of  the  scattering  operator, 
regardless  of  the  inversion  method  employed,  has  potential 
advantages  in  methods  of  collecting  and  analyzing  scattering 
data.  Previous  work  in  electrical  impedance  tomography  has 
employed  eigenfunction  decomposition  of  an  operator  asso¬ 
ciated  with  the  measurement  process  to  determine  optimal 
input  current  patterns  and  quantify  the  achievable  resolution 
of  imaging  systems.^’^  These  optimal  inputs  can  also  be  de- 
I  termined  by  iteratively  retransmitting  input  patterns  that  are 

proportional  to  the  measured  scattered  field.  This  approach  is 
essentially  an  analog  implementation  of  the  “power 

1^  method”  for  determining  the  eigenvectors  of  matrices.^’^ 

Likewise,  the  techniques  of  optical  and  acoustic  phase 
conjugation"^"^  and  the  analogous  process  of  time  reversal^’^ 


“turrent  affiliation:  Applied  Research  Laboratory,  The  Pennsylvania  State 
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can  be  understood  as  analog  methods  of  computing  the 
eigenfunctions  of  an  operator  associated  with  the  phase  con¬ 
jugation  or  time-reversal  process.  Simple  focusing  by  phase 
conjugation,  in  which  received  echoes  are  conjugated  or  time 
reversed  and  retransmitted,  is  equivalent  to  a  single  iteration 
of  the  power  method.  Further  iterations  of  this  procedure 
correspond  to  additional  steps  in  the  power  method,  and  thus 
converge  to  the  most  significant  eigenfunction  of  the  associ¬ 
ated  operator  at  a  rate  specified  by  the  ratio  of  the  two  largest 
eigenvalues.^  The  eigenfunctions  of  the  “time-reversal  op¬ 
erator,”  whether  obtained  by  iterative  time  reversal  or  by 
numerical  diagonalization,  have  been  previously  shown  to 
correspond  to  source  distributions  that  can  focus  incident 
energy  on  strong,  pointlike  scatterers.^^’^^ 

Eigensystem  analysis  has  historically  played  a  role  in 
the  theory  of  inverse  scattering  for  radially  symmetric 
objects. For  these  objects,  separation  of  variables  naturally 
leads  to  a  representation  of  the  scattering  operator  in  terms  of 
trigonometric  functions.  Since  these  eigenfunctions  are  the 
same  for  any  radial  scatterer,  the  inverse  scattering  problem 
could  be  reduced  to  the  problem  of  determining  the  unknown 
object  from  the  eigenvalues  of  the  scattering  operator. 
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However,  before  the  method  presented  here,  the  focus¬ 
ing  properties  of  eigenfunctions  have  not  been  exploited  for 
quantitative  reconstruction  of  inhomogeneous  media.  It  is 
stated  in  Ref.  9  that  the  concept  of  time  reversal  “cannot  be 
directly  compared  to  computed  tomography”  or  to  “tech¬ 
niques  that  generate  the  image  of  the  medium  through  signal 
analysis.”  Although  the  basic  principles  of  focusing  on  point 
targets  using  the  eigenfunctions  of  scattering  operators  have 
been  put  forth  in  Ref.  10,  these  principles  have  not  previ¬ 
ously  been  shown  to  apply  to  general  distributed  inhomoge¬ 
neities.  Furthermore,  no  general  imaging  method  has  hitherto 
been  based  on  these  principles. 

The  current  method  presents  a  solution  to  the  imaging 
problem  by  bringing  together  recent  results  in  the  theory  of 
focusing,  diffraction  tomography,  and  inverse  problems  to 
synthesize  a  unified  framework  for  quantitative  imaging  of 
inhomogeneous  media.  Application  of  the  method  shows  that 
focusing  on  distributed  inhomogeneities  can  be  achieved  us¬ 
ing  eigenfunctions  and  also  provides  a  technique  for  quanti¬ 
tative  imaging  of  discrete  and  distributed  inhomogeneities 
using  focusing  properties. 

This  method  has  several  advantages  over  current  inverse 
scattering  methods.  First,  the  eigenfunction  formulation  pro¬ 
vides  optimal  bases  for  reconstruction  of  unknown  media,  so 
that  inversions  are  performed  with  the  minimum  possible 
complexity.  Second,  the  method  is  applicable  to  any  scatter¬ 
ing  medium  for  which  the  total  acoustic  pressure  associated 
with  an  incident  plane  wave  can  be  estimated.  Inverse  scat¬ 
tering  in  inhomogeneous  background  media  as  well  as  itera¬ 
tive  nonlinear  inverse  scattering  can  therefore  be  directly 
implemented.  Third,  part  of  the  computation  necessary  for 
the  inverse  scattering  algorithm  can  be  performed  by  analog 
means  using  ideas  from  the  power  method. 

The  present  approach  also  provides  new  understanding 
about  existing  methods  of  focusing  and  imaging.  For  simple 
scattering  objects,  the  new  method  presented  here  reduces  to 
a  quantitative  specification  of  focusing  similar  to  that  ob¬ 
tained  by  iterative  phase  conjugation  or  time  reversal.  The 
eigenfunctions  of  scattering  operators  are  shown  not  only  to 
focus  on  pointlike  scatterers,  as  has  been  previously 
shown,^^’^^  but  also  to  concentrate  incident  energy  in  the 
vicinity  of  general,  distributed  inhomogeneities.  The  method 
also  improves  on  previous  approaches  to  focusing  using 
eigenfunctions  in  that  quantitative  images  of  medium  param¬ 
eters  are  obtained  simultaneously  with  optimal  incident- 
wave  distributions.  For  the  case  of  weakly  scattering  objects, 
the  method  reduces  to  a  simple  inversion  algorithm  that  is 
mathematically  equivalent  to  the  filtered  backpropagation 
algorithm,^^”^^  but  is  optimally  tailored  to  the  unknown  scat¬ 
tering  medium.  The  method  reduces  to  a  comparably  simple 
and  efficient  formula  for  the  case  of  weakly  nonlinear  in¬ 
verse  scattering. 

Analysis  given  in  Sec.  I  shows  that  eigenfunctions  of 
scattering  operators  are  equal  to  the  acoustic  fields  of  effec¬ 
tive  source  distributions  that  are  proportional  to  the  com¬ 
pressibility  contrast  of  the  scattering  object.  An  inverse  scat¬ 
tering  method  that  incorporates  products  of  retransmitted 
fields  of  eigenfunctions  is  presented.  The  general  method  is 
then  employed  to  derive  an  analytic  inversion  formula  valid 
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FIG.  1.  Scattering  configuration.  An  incident  plane  wave  traveling  in  the 
direction  a  is  scattered  by  an  inhomogeneity  and  the  scattered  field  is  mea¬ 
sured  in  the  direction  0. 


under  the  Bom  approximation  as  well  as  a  simple  nonlinear 
formula  valid  for  small  multiple-scattering  effects.  Numeri¬ 
cal  implementation  of  these  methods  is  presented  in  Sec.  II. 
Numerical  results  shown  in  Sec.  Ill  illustrate  focusing  on 
discrete  and  distributed  inhomogeneities  using  a  few  eigen¬ 
functions.  Also,  quantitative  inverse  scattering  results  are 
shown  both  within  the  context  of  a  homogeneous  back¬ 
ground  medium  and  an  inhomogeneous  background  medium. 

I,  THEORY 
A.  Background 

An  inverse  scattering  method  for  a  medium  of  variable 
sound  speed  is  derived.  For  simplicity  of  exposition,  the  deri¬ 
vation  is  given  for  the  canonical  two-dimensional  scattering 
configuration  sketched  in  Fig.  1.  However,  with  minor  modi¬ 
fications,  the  method  is  applicable  to  arbitrary  geometries 
and  dimensions. 

When  the  incident  pressure  is  a  plane  wave  of  unit  am¬ 
plitude  propagating  in  the  direction  a,  so  that 
—  the  corresponding  total  acoustic  pressure 

/?(x,a)  at  the  position  x  is  given  by  the  Lippman- Schwinger 
equation^^’^^ 

p(x,a)  =  e''“‘  j  Go(x-y,k)q(y)p(y,a)dy,  (1) 

where  Go(x— y)  is  the  Green’s  function  for  the  Helmholtz 
equation  in  a  homogeneous  medium.  In  an  unbounded  two- 
dimensional  medium,  Go(x-y)  is  given  by  the  Hankel  func¬ 
tion  (//4)Ho  \/:|x-y|).^^  The  angle  a  is  defined  as  the 
angle  corresponding  to  the  direction  unit  vector  a,  the  wave 
number  k  is  equal  to  lirffc^  where  Cq  is  the  wave  speed  of 
the  background  medium,  and  /  is  the  temporal  frequency  of 
the  incident  wave.  The  integral  appearing  in  Eq.  (1),  as  well 
as  subsequent  integrals  in  x  and  y,  are  understood  to  be 
taken  over  the  entire  plane  in  R^.  The  scattering  potential 
q  is  given  for  a  medium  of  variable  sound  speed  by 

The  quantity  within  parentheses  is  equal,  for  a  medium  of 
constant  density,  to  the  compressibility  contrast  ,  as  de¬ 
fined  in  Ref.  17.  The  scattering  potential  is  assumed  to  be 
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real-valued  and  to  be  short-range,  that  is,  the  potential  q 
decreases  at  large  distances  such  that 

|g(x)|^C(l-f|x|)-'-^  (3) 

where  |x|  is  the  magnitude  of  the  position  vector  x,  for  some 
S>0. 

At  a  measurement  radius  r  in  the  far  field  and  a  mea¬ 
surement  angle  6,  the  scattered  pressure,  Ps~P^Pi^  is  of 
the  form 

(4) 

where  A  is  the  far-field  pattern  of  the  scattered  pressure 

A{0,a)=  j  e~'^^'^q{x)p(x,a)dx.  (5) 

The  incident  pressure  may  be  more  generally  taken  as  a 
superposition  of  plane  waves  propagating  in  all  directions, 

/?j(x)=  J  da,  (6) 

The  far-field  pattern  of  the  corresponding  scattered  acoustic 
pressure  is  then 

J  A{6,a)f{a)da.  (7) 

Equation  (7)  defines  an  operator  A  that  maps  an  incident- 

wave  distribution  /( 0)  into  the  corresponding  far-field  scat¬ 
tered  pressure  Af{0).  The  operator  A  is  related  to  the  usual 
scattering  operator  S  (Ref.  19)  by 

S  =  I-~A,  (8) 

47r 

where  I  is  the  identity  operator. 

The  operator  A  is  compact^^  and  therefore  has  a  count¬ 
able  number  of  discrete  eigenvalues  with  zero  as  the  only 
possible  cluster  point.  In  practice,  only  a  finite  number  of 
eigenvalues  are  distinguishable  from  zero.  Since  the  poten¬ 
tial  q  is  real-valued,  the  scattering  operator  is  unitary,  so  that 
the  eigenvalues  of  A  lie  in  the  complex  plane  on  the  circle 
centered  at  -Am  and  passing  through  the  origin.  It  also 
follows  that  A  is  normal  (A*A=AA*,  where  A*  is  the  Her- 
mitian  transpose  of  A),  so  that  an  orthonormal  basis  {/,}  for 
L^[0,27r]  exists  consisting  of  eigenfunctions  of  A. 

Since  A  is  a  normal  operator,  the  Hermitian  transpose 
A*  satisfies  the  relation  A*fi  =  K*fi,  where  fi  is  an  eigen¬ 
function  of  A  and  X*  is  the  complex  conjugate  of  X^.  The 
eigenfunctions  of  A  therefore  also  satisfy  the  equation 

(9) 

Thus  the  functions/^-  also  constitute  a  basis  of  eigenfunctions 
for  A  *A  and  the  corresponding  eigenvalues  are  the  squared 
magnitudes  of  the  eigenvalues  of  A.  The  operator  A*A  is 
essentially  a  far-field  analog  of  the  “time-reversal  operator” 
as  defined  in  Ref.  10. 


B.  Focusing  properties 

The  focusing  properties  of  A  are  seen  by  considering  the 
ratio  of  the  scattered  amplitude  to  the  incident  amplitude. 
Since  A  is  normal,  the  magnitude  of  its  largest  eigenvalue  is 
equal  to  the  largest  possible  value  of  this  ratio  for  any  non¬ 
zero  /: 


lXi|  =  sup 


II/(^)IIl2  , 


(10) 


where  sup(‘)  denotes  the  least  upper  bound  and  li/(-)llL2 
denotes  the  root-mean-square  magnitude  of  a  square- 
integrable  function.  Thus  the  eigenfunction  associated  with 
the  largest  eigenvalue  of  A  specifies  an  incident-wave  distri¬ 
bution  that  maximizes  the  energy  scattered  to  the  far  field. 
Other  eigenfunctions  also  focus  energy  on  inhomogeneities 
with  an  efficiency  that  is  quantified  by  the  associated  eigen¬ 
values. 

The  focusing  property  of  eigenfunctions  of  A  is  further 
illustrated  by  introducing  the  acoustic  fields  of  incident-wave 
distributions  specified  by  the  eigenfunctions.  One  may  define 
retransmitted  fields  of  an  incident- wave  distribution /(a)  as 


£(x)=  da, 

F(x)=J  f(a)p(x,a)da, 


(11) 


where  E(x)  is  the  retransmitted  field  associated  with  the 
incident-wave  distribution  in  a  homogeneous  medium  and 
F(x)  is  the  retransmitted  field  in  a  medium  containing  the 
inhomogeneity  ^(x). 

For  incident-wave  patterns  corresponding  to  eigenfunc¬ 
tions  that  have  nonzero  eigenvalues,  the  retransmitted  fields 
of  Eqs.  (11)  can  be  written  using  Eqs.  (5)  and  (7)  in  the  form 


Ei(x)=  Jo(klx-yl)Fi(y)q(y)dy, 

The  brackets  in  Eq.  (12)  denote  the  inner  product 


(12) 


u{d)v*{0)dd. 


(13) 


while  the  inner  product  appearing  in  the  expression  for  E, 
was  evaluated  using  the  identity 

(14) 

iTT  J 

known  as  Parseval’s  integral.^®  The  retransmitted  fields  of 
Eq.  (11)  are  thus  seen  to  be  equivalent  to  a  weighted  convo¬ 
lution  of  the  unknown  scattering  potential  with  inner  prod¬ 
ucts  of  acoustic  fields. 

When  the  scattering  potential  ^(x)  is  concentrated  in  a 
finite  number  of  pointlike  scatterers,  each  very  small  com¬ 
pared  to  a  wavelength,  Eq.  (12)  reduces  to  an  expression  of 
diffraction-limited  focusing  on  each  point  scatterer.  That  is, 
for  a  scattering  medium  defined  by 
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M 

?(x)=S  /U.;5(x-Xj),  (15) 

the  retransmitted  field  ^/(x)  is 

£,(x)=^2  £,(xy)yo(^|x-Xj|)/i;,  (16) 

j 

SO  that  in  this  case,  the  retransmitted  field  ^/(x)  is  equal  to  a 
weighted  sum  of  Bessel  functions,  each  centered  at  the  loca¬ 
tion  of  one  of  the  point  scatterers.  These  Bessel  functions 
correspond  to  a  group  of  diffraction-limited  main  lobes,  cen¬ 
tered  at  each  scatterer  position  Xj,  with  corresponding  Bessel 
sidelobes  that  combine  coherently.  Thus  each  retransmitted 
field  El  focuses  to  some  extent  on  all  of  the  individual  point 
scatterers. 

The  close  relationship  between  the  retransmitted  fields 
of  eigenfunctions  and  the  unknown  scattering  potential,  as 
seen  in  Eq.  (12),  is  an  expression  of  the  focusing  property  of 
eigenfunctions.  That  is,  since  eigenfunctions  of  A  correspond 
to  incident-wave  patterns  that  concentrate  energy  within  the 
support  of  the  scattering  potential,  they  can  be  said  to  focus 
on  general  distributed  inhomogeneities  as  well  as  pointlike 
scatterers.  This  idea  is  illustrated  numerically  later  in  this 
paper. 


C.  Inverse  scattering  method 

Because  of  the  focusing  properties  outlined  above,  re¬ 
transmitted  fields  of  eigenfunctions  are  a  useful  starting  point 
for  inverse  scattering  reconstructions.  A  general  inverse  scat¬ 
tering  method  incorporating  these  ideas  is  outlined  below. 

The  starting  point  for  this  method  is  an  expression  of  the 
inverse  scattering  problem  in  terms  of  the  operator  A  of 
Eq.  (7)  and  the  corresponding  retransmitted  fields  defined  in 
Eq.  (11): 

{Afi ,fj)=-Sij\i=  j  Fi(x)EJ(x)^(x)dx, 

ij=  1,2,...  .  (17) 

The  problem  can  be  regularized  by  seeking  the  solution  that 
minimizes  the  weighted  norm 

\q{x)\^'W{x)dx  (18) 

with  lV(x)  an  appropriate  weight.  For  the  analysis  given  be- 
low,  this  weight  is  defined  as  iy(x)  =  ( 1  +  Ix])*^,  <5>0.  For  the 
explicit  computations  given  later,  other  choices  of  W(x)  are 
more  natural. 

A  solution  to  the  minimization  problem  is  obtained  us¬ 
ing  the  method  of  Lagrange  multipliers,  analogous  to  the 
approach  used  in  Ref.  21  for  a  linearized  electric  impedance 
tomography  problem.  At  a  minimum,  the  (infinite¬ 
dimensional)  gradient  of  II  ^11^  is  a  linear  combination  of  the 
gradients  of  the  constraints  in  Eq.  (17).  The  latter  can  be 
calculated  using  the  two-potential  formula^^ 


Aq^{e,a)-A^^{e,a)  = 


/?i(x,a) 


X(^i(x)-^2(x))P2(x,6>+7r)^/x, 

(19) 

where  px ,  A^^,  and  pi  are  the  scattering  operators  and 
the  total  acoustic  pressures  for  the  inhomogeneous  media 
defined  by  ^i(x)  and  qiix),  respectively.  Equation  (19) 
yields  the  derivative 


lim 

e--+0 


Aq  +  e2i{0,a)-A^{e,a) 


=  J  p{x,a)q{x)p{x,d+  'w)dx, 

(20) 


while  the  infinite-dimensional  gradient  of  ll^ll^  is  found  from 


lim 


\\q  +  €q\\\- 


■WqWl 


=  2  j  q(x)q(x)W(x)dx.  (21) 


The  result  follows  that  if  the  potential  qui^)  solves  the 
regularized  inverse  scattering  problem  [minimization  of  the 
weighted  norm  from  Eq.  (18)  under  the  constraint  of  Eq. 
(17)],  qM  must  be  of  the  form 

1^77^2)  2  (22) 

W(Xj  /  m 


where  F^(x),  the  complex  conjugate  of  the  retransmitted 
field  corresponding  to  an  incoming  condition  at  infinity,  is 
defined  as 


f*(a)p(x,a+  7r)da, 


(23) 


and  the  coefficients  Qi^  are  the  Lagrange  multipliers.  If  the 
above  gradients  are  taken  with  respect  to  the  real  and  imagi¬ 
nary  parts  of  a  complex  potential,  Eq.  (22)  as  stated  is  also 
found  to  be  valid  when  the  potential  is  complex.  In  some 
of  the  simplifying  approximations  made  below,  Eq.  (22)  will 
yield  a  complex  potential  q^^  even  when  the  data  are  as¬ 
sumed  to  come  from  a  unitary  scattering  operator  associated 
with  the  real  potential  q. 

By  substituting  Eq.  (22)  into  Eq.  (17),  the  inverse  prob¬ 
lem  is  reduced  to  the  problem  of  finding  the  coefficients 
Qi^  from  the  nonlinear  system 


2  f 

^  I  m  J 


Fi(x)E]ix)Fiix)F*Jx) 

W{x) 


dx 


Q 


Im  ’ 


£7  =  1,2,...,  (24) 

where  the  dependence  of  the  fields  F  and  F*  on  the  scatter¬ 
ing  potential  q  is  implicit. 

In  general,  the  scattering  potential  q{x),  and  therefore 
the  total  pressure  field  p(x,a),  are  unknown  in  inverse  scat¬ 
tering  problems.  The  function  /?(x,a)  that  implicitly  appears 
in  Eq.  (24)  may  therefore  be  replaced  by  the  best  available 
estimate  for  the  total  pressure.  Equation  (24)  can  then  be 
solved  for  the  coefficients  Qi^  by  standard  numerical  tech¬ 
niques  for  solution  of  linear  systems. 

The  number  of  terms  N  can  be  chosen  arbitrarily;  how¬ 
ever,  increasing  N  beyond  the  number  of  nonzero  eigenval- 
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ues  of  A  is  of  limited  benefit  in  reconstructions.  For  simple 
scattering  objects,  q  can  be  represented  by  expansions  em¬ 
ploying  small  values  of  N.  For  instance,  for  an  inhomogene¬ 
ity  consisting  of  finitely  many  point  scatterers,  N  comparable 
to  the  number  of  scatterers  is  sufficient. 

The  above  method  simplifies  further  in  the  case  of  a 
weakly  scattering  medium,  for  which  the  total  pressure  p  can 
be  approximated  by  the  incident  pressure.  In  this  case,  taking 
the  weight  lT(x)  =  l,  the  coefficients  can  be  evaluated 
analytically.  From  Eq.  (22),  under  the  Bom  approximation, 
the  scattering  potential  takes  the  form 

?5(x)  =  S  E  Qt„Ei{x)E*Jx).  (25) 

l  m 

Substituting  Eq.  (25)  into  Eq.  (5)  gives  the  equation 

A(0,ff)  =  2  S  Qiml  e~''^^'^Ei{x)E*Jx)p{x,a)dx. 

I  m  J 

(26) 

Replacement  of  p(x,a)  in  Eq.  (26)  by  the  incident  plane 
wave  use  of  Eq.  (11),  and  integration  in  x  over 

yields 

A{0,a)=-p-^J,Q,„\  S{e-a-e'  +  a') 

K  I  m  J  —  TT  J  —  TT 

'Xfi{e^)fl{a^)d0^da^  (27) 

for  a  not  equal  to  0  or  tt.  The  double  integral  in  Eq.  (27) 
can  be  evaluated  using  the  change  of  variables 

x|  =  cos  — cos  a\  A:2  =  sin  — sin  a\  (28) 

which  is  one-to-one  when  restricted  to  the  regions  a'<6' 
and  a^>0' .  Evaluation  of  the  integral  yields 

(2,7r)^ 

|sin(0-a)|A(0,a)=— — E  E  QimifiWU^) 
k  I  tn 

+//(«+ 77)/*  (^+7r)).  (29) 

Equation  (29)  can  be  solved  for  the  coefficients  Qi^  using 
the  fact  that  the  eigenfunctions //( 9)  are  orthonormal  as  well 
as  the  reciprocity  identity^^ 

A  ( 77,  a  +  77)  =  A  ( a,  ^) ,  (30) 

The  solution  is 

=  oi)\Aie,a)fJ{e)fJa)da  dd. 

(31) 

Equations  (25)  and  (31)  specify  a  solution  to  the 
linearized  inverse  problem.  This  solution  is,  in  general,  com¬ 
plex,  even  when  the  tme  potential  q  is  purely  real.  A  physi¬ 
cal  way  to  understand  why  the  Bom  approximation  yields  a 
complex  scattering  potential  for  a  lossless  medium  is  to  rec¬ 
ognize  that  this  approximation  neglects  multiple  scattering 
and  thus,  the  resulting  output  energy  differs  from  the  input 
energy.  The  corresponding  scattering  operator  is  then  no 
longer  unitary,  and  is  only  physically  realizable  by  a  poten¬ 
tial  with  a  nonzero  imaginary  part.  For  weak  scattering,  the 
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energy  discrepancy  is  small  and  so  is  the  imaginary  part  of 
the  potential. 

The  analytic  solution  of  Eqs.  (25)  and  (31)  is  equivalent 
to  the  well-known  filtered  backpropagation  formula^^"^^  and 
has  the  advantage  of  computational  simplicity,  as  discussed 
later  in  this  paper.  Equivalence  between  the  two  formulas  is 
shown  by  formulating  an  expansion  of  viewed  as  a 

function  of  a,  in  terms  of  the  orthonormal  basis  {/m(<^)}-  lo 
view  of  Eqs.  (11),  this  expansion  yields  the  identity 

£:(x)/.(a).  (32) 

m 

Substituting  Eq.  (31)  in  Eq.  (25)  and  using  Eq.  (32)  as  well 
as  its  conjugate  gives 

^B(x)=^J|lsin(a-^)lA(0,a)e'**  ('’-“Va  dd,  (33) 

which  is  the  standard  filtered  backpropagation  formula. 
Equation  (33)  yields  the  low-pass-filtered  version  of  the  tme 
potential  q  if  multiple  scattering  effects  are  negligible.  The 
correct  nonlinear  generalization  of  the  linearized  low-pass 
filtered  solution  ^5  is  the  minimal  (or  weighted  L^)  so¬ 
lution  qj^ ,  which  is  of  a  form  specified  by  Eq.  (22). 

The  inverse  scattering  method  developed  above  can  also 
be  used  with  any  orthonormal  set  of  basis  functions  for 
L^[0,27r].  For  instance,  reconstmctions  can  be  performed 
using  eigenfunctions  of  A  for  axisymmetric  objects  rather 
than  using  the  eigenfunctions  associated  with  the  measured 
A .  In  this  case,  the  eigenfunctions  take  the  form 

m  =  0,±l,±2,...  .  (34) 

V277 

The  retransmitted  fields  can  be  analytically  evaluated  to 
be 

EJr,<f,)=^j2^re‘'"‘l’JJkr),  (35) 

and  the  coefficients  Qi^  for  the  low-pass-filtered  reconstmc- 
tion  of  q  are  given  by 

Qin,-^ / / \Me-a)\Aie,a)e-“^e‘"‘“da  dS. 

(36) 

While  the  retransmitted  fields  specified  by  Eq.  (35)  are  not 
ideally  matched  to  nonaxisymmetric  scattering  media,  they 
can  be  analytically  evaluated  and  stored  for  use  in  fast  re¬ 
constmctions.  Since  these  retransmitted  fields  are  also  unaf¬ 
fected  by  uncertainties  in  scattering  measurements,  they  are 
suitable  for  reconstmctions  from  data  cormpted  by  noise. 

Finally,  use  of  the  eigenfunction  method  beyond  linear 
inversion  is  demonstrated  by  considering  the  case  where  the 
inhomogeneous-medium  retransmitted  fields  F  can  be  esti¬ 
mated  from  a  first  approximation  to  the  scattering  potential 
q.  One  approach  in  this  case  is  to  solve  the  full  system  of 
equations  defined  by  Eq.  (24);  however,  a  more  numerically 
efficient  correction  to  the  Bom  approximation  can  be  ob¬ 
tained  by  invoking  the  localized  nonlinear  approximation  in¬ 
troduced  in  Ref,  22  for  electromagnetic  scattering.  This  ap¬ 
proximation  follows  from  writing  the  Lippman-Schwinger 
equation  [Eq.  (1)]  in  the  form 
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/j(x,a)  =  r(x)|e'^“''’‘-  j  (p(y,a)-p(x,a)) 

X<l(y)Go(^-y)dyj,  (37) 

where  the  quantity  r(x),  called  the  depolarization  tensor  in 
electromagnetic  scattering, is  defined  by 

r(x)  =  |n- J  ^(y)Go(x-y)Jy|  .  (38) 

The  second  term  in  Eq.  (37)  is  presumed  to  be  small  because 
the  singularity  of  the  Green’s  function  is  cancelled  by  the 
difference  term  appearing  in  the  integrand.  Thus  the  total 
pressure  may  be  approximated  by  the  formula 

p(x,a)^r(x)e^^^'\  (39) 

The  form  for  the  scattering  potential  given  by  Eq.  (22) 
then  becomes 

^Ar(x)-^S  E  Q,„Mx)E:(x).  (40) 

Substituting  this  form  into  Eq.  (5)  and  using  Eq.  (39) 
gives 

I  m  J 

xr{x)^e‘‘^“-’‘ dx.  (41) 

An  approximate  nonlinear  formula  for  the  scattering  poten¬ 
tial  q  can  be  obtained  by  taking  lT(x)  =  r(x)^.  Equation  (41) 
then  yields  the  coefficients  from  Eq.  (31)  and  the  result¬ 
ing  solution  for  the  scattering  potential  is 

_  ^  ^  ^  E,{x)E*Jx) 

2  Qlm  p/  \  •  (42) 

I  m  1  (Xj 

The  solution  of  Eq.  (42)  is  simplified  by  making  the 
further  approximation 

f^-2-r(x).  (43) 

which  is  valid  for  small  scattering  potentials.  This  substitu¬ 
tion  results  in 

?a/(x)-E  E  QiME,{x)-F,ix))Et,{x).  (44) 

/  m 

This  nonlinear  equation  for  the  potential  can  be  approxi¬ 
mately  solved  by  using  a  form  of  the  retransmitted  field 
F/(x)  corresponding  to  the  low-pass-filtered  potential  qs  or 
to  another  estimate  of  the  scattering  potential. 

II.  COMPUTATIONAL  METHODS 

The  focusing  and  imaging  methods  outlined  in  Sec.  I 
were  implemented  using  numerically  computed  scattered 
fields  of  inhomogeneous  objects.  Scattering  operators  were 
calculated  using  a  method  due  to  Kirsch  and  Monk,^^  in 
which  an  inner  solution  of  the  Helmholtz  equation  for  a  me¬ 
dium  of  variable  sound  speed  is  matched  to  an  outer  solution 
of  integral  equations  that  implicitly  satisfy  the  Sommerfeld 
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radiation  condition.  The  inner  solution  is  obtained  using  a 
finite-element  method,  while  the  outer  integral  equations  are 
solved  using  Nystrom’s  method.^"^ 

Scattering  data  were  calculated  numerically  for  a  num¬ 
ber  of  incident  plane  waves  evenly  distributed  over  M  angles 
between  0  and  2  tt.  For  each  incident-wave  angle,  the  scat¬ 
tered  field  was  computed  at  M  far-field  receiver  angles  be¬ 
tween  0  and  277,  so  that  the  angular  sampling  rate  was 
Mfilrr)  samples  per  radian.  The  number  of  receiver  angles 
M  should  be  chosen  such  that  the  scattered  field  has  no  sig¬ 
nificant  frequency  components  above  the  Nyquist  frequency 
of  M/(47r)  samples  per  radian.  This  computation  yields  a 
discrete  representation  of  the  scattering  operator  A  as  an 
MXM  matrix,  A ^ . 

The  eigenfunctions  of  A  and  their  associated  eigenval¬ 
ues  were  estimated  numerically  by  direct  computation  of  the 
eigenvectors  and  eigenvalues  of  .  Retransmitted  fields  of 
eigenfunctions  were  evaluated  numerically  by  numerical  in¬ 
tegration  of  Eq.  (11).  Images  of  inhomogeneous  objects  were 
then  obtained  using  a  straightforward  numerical  implemen¬ 
tation  of  Eqs.  (25)  and  (31).  The  integrals  appearing  in  these 
equations  were  evaluated  using  corresponding  discrete  sum¬ 
mations  of  the  components  of  A  ^  and  its  eigenvectors.  For 
comparison,  standard  diffraction  tomography  inversions 
were  also  performed  by  numerical  integration  of  Eq.  (33). 

Stability  of  the  eigenfunction  imaging  method  was 
tested  by  inversion  of  noisy  data  obtained  by  adding  numeri¬ 
cally  generated  Gaussian  white  noise  to  the  scattering  matrix 
A^ .  The  rms  amplitude  of  the  noise  was  specified  as  a  frac¬ 
tion  of  the  rms  value  of  A  .  Thus,  for  instance,  a  signal-to- 
noise  ratio  of  6  dB  was  obtained  by  adding  noise  with  an  rms 
amplitude  one-half  the  rms  value  of  A  • 

Inversions  were  also  performed  using  the  basis  of  eigen¬ 
functions  corresponding  to  axisymmetric  inhomogeneities. 
In  this  case,  the  formula  of  Eq.  (25)  was  implemented  nu¬ 
merically  using  the  trigonometric  basis  functions  defined  in 
Eq.  (34),  the  retransmitted  fields  given  in  closed  form  in  Eq. 
(35),  and  the  coefficients  defined  in  Eq.  (36). 

Nonlinear  eigenfunction  images  were  obtained  using  the 
analytic  formula  of  Eq.  (44)  with  the  total  pressure  p(x,a) 
approximated  by  the  total  pressure  for  a  medium  including  a 
cylinder  of  specified  radius  and  compressibility  contrast. 
This  computation  employed  an  exact  solution  for  the  scatter¬ 
ing  of  a  plane  wave  by  a  cylinder. 

III.  NUMERICAL  RESULTS 

Focusing  of  eigenfunctions  on  a  distributed  scattering 
object  is  illustrated  in  Fig.  2.  Here,  the  magnitudes  of  the 
retransmitted  fields  £:i(x)  and  E2(x)  are  shown  for  an  inho¬ 
mogeneity  consisting  of  a  weakly  scattering  triangle 
(7^=0.01)  approximately  two  wavelengths  in  height.  The 
triangle  is  shown  in  outline  together  with  the  retransmitted 
fields.  The  corresponding  scattering  operator,  calculated  us¬ 
ing  the  finite-element/Nystrom  method  described  above,  was 
represented  by  a  matrix  of  size  128X  128.  The  retransmitted 
fields  show  that  the  significant  eigenfunctions  of  A  specify 
incident-wave  patterns  that  concentrate  energy  within  the 
support  of  the  inhomogeneity.  Notable  is  that  this  focused 
energy  is  distributed  throughout  the  support  of  the  triangle. 
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Implementation  of  the  eigenfunction  method  in  focusing 
on  pointlike  scatterers  is  illustrated  in  Figs.  3  and  4,  These 
figures,  obtained  using  the  linearized  eigenfunction  method, 
show  not  only  diffraction-limited  focusing  but  also  quantita¬ 
tive  reconstructions  of  the  associated  scattering  potentials. 
Figure  3  shows  images  made  from  the  scattered  field  of  two 
pointlike  scatterers  at  locations  (  —  0.5,0)  and  (0,— 0.2),  each 
of  radius  0.01  and  compressibility  variation  -0.9.  The  nu¬ 
merically  computed  scattered  field  was  sampled  at  128 
equally  spaced  angles  for  each  of  128  incident- wave  angles, 
so  that  the  operator  A  was  represented  by  a  128  X  128  matrix. 
The  wave  number  used  was  10,  so  that  the  scatterers  were 
separated  by  approximately  one  wavelength.  Since,  in  this 
case,  two  eigenvalues  of  A  were  much  larger  than  the  re¬ 
maining  eigenvalues,  the  basic  reconstruction  required  only 
the  use  of  two  retransmitted  fields.  This  result  illustrates  that, 
for  an  inhomogeneity  consisting  of  finitely  many  pointlike 
scatterers,  the  present  inverse  scattering  method  provides  an 
accurate  reconstruction  with  diffraction-limited  point  resolu¬ 
tion  using  a  corresponding  number  of  eigenfunctions. 

A  stability  test  of  the  eigenfunction  method  is  illustrated 
in  Fig.  4.  This  image  shows  a  reconstruction  of  the  two 
pointlike  scatterers  of  Fig.  3  using  the  same  scattering  data 
with  added  Gaussian  white  noise  for  a  signal-to-noise  ratio 
of  3  dB.  The  method  of  reconstruction  was  identical  to  that 
used  for  Fig.  3.  The  reconstruction  shown  is  almost  indistin¬ 
guishable  from  the  noiseless  reconstruction,  indicating  the 
stability  of  the  eigenfunction  imaging  method. 

Linear  eigenfunction  images  of  the  triangular  inhomoge¬ 
neity  of  Fig.  2  are  presented  in  Fig.  5.  These  images  were 
constructed  using  the  same  scattering  data  as  that  used  for 
Fig.  2.  The  first  image,  obtained  using  five  retransmitted 
fields,  shows  that  strong  focusing  is  achieved  using  only  a 
few  eigenfunctions  of  A .  The  entire  inhomogeneity  is  well- 
insonified  and  little  incident  energy  is  transmitted  outside  the 
support  of  the  inhomogeneity.  The  second  image,  obtained 
using  15  eigenfunctions,  shows  that  the  eigenfunction 
method  rapidly  converges  to  the  ideal  low-pass-filtered  solu¬ 
tion  for  the  scattering  potential.  Notable  is  that  the  eigen¬ 
function  method  using  15  eigenfunctions  required  69.1  s  of 
CPU  time  on  a  Sun  SPARCstation  10,  while  an  analogous 
image  obtained  using  the  diffraction  tomography  formula  of 
Eq.  (33),  with  the  integrals  evaluated  in  an  analogous  man¬ 
ner,  required  3014.3  s. 

Eigenfunction  reconstructions  of  a  test  phantom,  shown 
in  Figs.  6-8,  illustrate  application  of  the  eigenfunction  im¬ 
aging  method  to  a  larger-scale  imaging  problem.  The  phan¬ 
tom,  also  represented  in  Fig.  1,  is  a  cylinder  of  compressibil¬ 
ity  contrast  0.01  and  diameter  of  5  mm.  Internal  objects 
include  a  water-filled  (cystic)  region  of  diameter  1  mm,  a 
wire  of  diameter  0.1  mm  and  compressibility  contrast 
—  0.5,  and  an  internal  cylinder  of  diameter  1  mm  and  com¬ 
pressibility  contrast  -0.01.  Scattered  fields  were  calculated 
using  the  methods  described  above,  with  the  operator  A  dis¬ 
cretized  as  a  matrix  of  256X256  points.  The  first  image 
shown  in  Fig.  6,  obtained  using  the  single  wave  number 
A:==10  has  high  resolution  but  contains  ringing  (Gibbs  phe¬ 
nomenon)  artifacts  and  loss  of  contrast  in  the  cystic  region. 
These  artifacts  are  removed  by  compounding  of  images  ob- 
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tained  using  five  linearly  spaced  wave  numbers,  8^/:^  12, 
so  that  the  dimensionless  parameter  ka  varied  between  20 
and  30.  The  five-frequency  image,  shown  in  the  second  panel 
of  Fig.  6,  also  shows  increased  point  and  contrast  resolution 
compared  to  the  single-frequency  image.  Both  images  shown 
in  Fig.  6  were  obtained  using  the  linearized  eigenfunction 
method  described  above,  with  64  eigenfunctions  of  A  for 
A:=8,  9,  and  10,  68  eigenfunctions  for  A:  =  11,  and  72  eigen¬ 
functions  for  A;=  12. 

Reconstructions  of  the  test  phantom  obtained  from  noisy 
data  are  shown  in  Fig.  7.  Gaussian  white  noise  was  added  to 
the  scattering  data  employed  in  Fig.  6,  so  that  the  signal-to- 
noise  ratio  was  6  dB  at  each  of  the  frequencies  employed. 
The  reconstructions  employed  the  formula  of  Eq.  (25)  and 
coefficients  obtained  from  Eqs.  (34)-(36).  The  numbers  of 
basis  functions  employed  were  equal  to  the  number  of  eigen¬ 
functions  employed  in  Fig.  6.  These  results  indicate  the  sta¬ 
bility  of  the  method  for  large  objects  with  scattering  data 
severely  degraded  by  noise. 

Nonlinear  reconstructions  of  the  same  test  phantom,  ob¬ 
tained  using  Eq.  (44),  are  presented  in  Fig.  8  together  with 
linear  reconstructions.  In  the  nonlinear  reconstructions,  the 
retransmitted  fields  F/(x)  were  estimated  using  pressure 
fields  associated  with  a  cylinder  of  diameter  5  mm  and  com¬ 
pressibility  contrast  0.01.^^  The  scattering  data  employed 
was  identical  to  that  used  in  Fig.  6(b),  with  five  linearly 
spaced  wave  numbers  such  that  20^ka^30.  The  number  of 
eigenfunctions  employed  in  each  image  were  also  the  same 
as  those  used  for  the  images  in  Fig.  6.  The  first  panel  shows 
the  real  part  of  the  nonlinear  reconstruction,  taken  along  the 
line  y  =  0,  together  with  the  real  part  of  the  analogous  linear 
reconstruction  from  Fig.  6(b).  The  nonlinear  reconstruction 
shows  improved  resolution  over  the  linear  reconstruction  by 
increased  height  of  the  peak  associated  with  the  internal 
wire.  The  second  panel  shows  the  imaginary  part  of  the  non¬ 
linear  reconstruction  with  the  corresponding  linear  recon¬ 
struction  from  Fig.  6(b).  Here,  the  inaccuracy  of  the  Bom 
approximation  results  in  a  significant  imaginary  part  for  the 
linear  reconstruction,  while  the  tme  potential  is  purely  real. 
The  nonlinear  inversion  shows  improved  quantitative  accu¬ 
racy  over  the  linear  inversion  by  reduction  of  the  recon¬ 
structed  imaginary  part. 

IV.  DISCUSSION 

Our  method  has  shown  that  eigenfunctions  of  the  scat¬ 
tering  operator  can  be  employed  to  focus  on  distributed  in¬ 
homogeneities  as  well  as  pointlike  scatterers.  However,  the 
focusing  on  distributed  inhomogeneities  occurs  in  a  different 
manner  from  focusing  on  pointlike  scatterers.  That  is,  the 
incident  energy  is  not  maximized  at  a  single  point  within  the 
medium.  Instead,  when  combined  according  to  the  present 
reconstmction  method,  retransmitted  eigenfunctions  specify 
incident-wave  distributions  that  distribute  energy  throughout 
the  inhomogeneous  region.  This  type  of  focusing,  which  re¬ 
sults  from  the  eigenfunction  property  of  maximizing  the 
scattered  energy,  is  clearly  connected  to  imaging  of  the  me¬ 
dium  by  inverse  scattering. 

The  quantitative  inverse  scattering  method  presented  in 
this  paper  can  considerably  simplify  imaging  computations. 
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FIG.  2.  Focusing  on  a  distributed  inhomogeneity.  Magnitudes  of  the  retransmitted  fields  of  the  two  most  significant  eigenfunctions  are  shown  on  a  linear  gray 
scale  with  black  indicating  zero  and  white  indicating  maximum  amplitude.  The  scattering  object  is  a  uniform  triangle,  compressibility  contrast  0.01,  within  the 
sketched  boundaries,  (a)  (b)  £’2(x). 


The  computational  complexity  of  the  current  method  de¬ 
pends  mainly  on  the  number  of  significant  eigenfunctions, 
which  in  turn  depends  only  on  the  complexity  of  the  scatter¬ 
ing  medium.  Furthermore,  the  basis  for  expansion  of  the  un¬ 
known  medium  is  determined  directly  from  the  scattering 
data.  Since  the  basis  functions  are  closely  related  to  the  un¬ 
known  medium,  reconstructions  performed  using  this  basis 
employ  a  minimal  amount  of  unnecessary  information.  This 
property  gives  the  present  inverse  scattering  method  advan¬ 
tages  over  other  methods  in  which  a  fixed  basis  is  used  to 
expand  the  unknown  medium.^"^"^^  These  advantages  are 


most  apparent  for  inhomogeneities  a  small  number  of  wave¬ 
lengths  in  size. 

The  present  inverse  scattering  method  also  has  the  ad¬ 
vantage  of  applicability  to  any  medium  for  which  the  back¬ 
ground  pressure  field  can  be  estimated.  Use  of  background 
pressure  estimates  can  greatly  improve  accuracy  over  recon¬ 
structions  based  on  simpler  approximations.  For  instance, 
Born  inversion  can  yield  a  spurious  reconstructed  imaginary 
part  even  when  the  true  potential  is  real-valued;  use  of  an 
estimated  background  pressure  field  can  greatly  reduce  this 
error,  as  seen  in  Fig.  8.  The  inverse  scattering  method  pre- 


FIG.  3.  Eigenfunction  image  of  two  pointlike  scatterers,  compressibility 
contrast  -0.9,  separated  by  approximately  one  wavelength.  The  image  was 
obtained  using  retransmitted  fields  of  the  two  most  significant  eigenfunc¬ 
tions. 


FIG,  4.  Effect  of  noise  on  eigenfunction  reconstruction.  The  object  of  Fig.  3 
was  reconstructed  from  two  eigenfunctions  of  synthetically  noised  scattering 
data  with  a  signal-to-noise  ratio  of  3  dB. 
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FIG.  5.  Eigenfunction  images  of  a  triangle  about  three  wavelengths  in  height  having  compressibility  contrast  0.01.  (a)  Real  part  of  an  inversion  obtained  using 
retransmitted  fields  of  five  eigenfunctions,  (b)  Analogous  image  obtained  using  15  eigenfunctions. 


sented  here  is  also  extensible  to  any  background  medium  for 
which  a  pressure  field  can  be  estimated,  including  moving 
fluids,  layered  or  stratified  media,  and  enclosed  or  otherwise 
bounded  regions. 

The  imaging,  focusing,  and  inverse  scattering  methods 
presented  here  also  intrinsically  take  advantage  of  any  poten¬ 
tial  increase  in  resolution  that  is  associated  with  multiple 
scattering  or  other  higher-order  effects,  as  long  as  these  ef¬ 
fects  are  taken  into  account  in  the  estimated  pressure  field. 
This  increase  in  resolution  is  possible  because  the  retransmit¬ 
ted  fields  of  eigenfunctions  may  have  desirable  qualities, 
such  as  higher  spatial-frequency  components,  that  are  asso¬ 
ciated  with  the  presence  of  an  inhomogeneous  background. 
Such  improvements  in  resolution  have  been  shown  experi¬ 
mentally  for  time-reversal  focusing  in  a  multiply  scattering 
medium.^^ 


The  most  significant  eigenfunctions  of  A  can  be  deter¬ 
mined  experimentally  through  iterative  retransmission  of  re¬ 
ceived  scattered  fields  in  a  manner  similar  to  that  performed 
in  Refs.  1  and  10.  This  implementation  of  the  power  method^ 
allows  computation  of  the  most  significant  eigenfunctions  of 
A  by  analog  means,  which  may  be  preferable  to  numerical 
computation  for  very  large  scattering  objects.  These  eigen¬ 
functions  are  useful  as  optimal  incident-wave  patterns  for 
inverse  scattering  experiments. 

The  inverse  scattering  method  presented  here  includes 
the  assumption  that  the  scattering  potential  is  purely  real, 
that  is,  the  inhomogeneous  medium  is  assumed  to  have  zero 
absorption.  Eigensystem  analysis  of  the  scattering  operator 
A  is  more  complicated  in  the  presence  of  absorption. 
However,  the  methods  introduced  here  are  expected  to  still 


FIG.  6.  Eigenfunction  images  of  a  test  object  having  background  compressibility  contrast  0.01,  a  cystic  (water-filled)  region,  a  pointlike  scatterer,  and  an 
internal  cylinder.  The  images  are  shown  on  a  logarithmic  scale  with  a  40  dB  dynamic  range,  (a)  Real  part  of  inversion  obtained  for  a  wave  number  such  that 
ka-25.  (b)  Analogous  image  obtained  using  five  wavenumbers  such  that  20<ka<30. 
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FIG.  7.  Images  of  test  object  obtained  using  noised  data  with  6  dB  signal -to-noise  ratio.  The  images  are  shown  on  a  logarithmic  scale  with  a  40  dB  dynamic 
range,  (a)  Real  part  of  inversion  obtained  for  a  wave  number  such  that  ka  =  25.  (b)  Analogous  image  obtained  using  five  wave  numbers  such  that  20<ka 
<30. 


be  useful  for  media  such  as  biological  tissue  when  absorption 
is  finite  but  small. 

A  disadvantage  of  the  inverse  scattering  method  as  cur¬ 
rently  implemented  is  that  nonlinear  inversion  requires  an 
accurate  specification  of  the  background  acoustic  field.  This 
disadvantage  is  not  unique  to  the  eigenfunction  method,  but 
is  a  common  feature  of  most  current  inversion  methods.  Re¬ 
cent  theoretically  exact  methods,  while  not  limited  in  this 
manner,^ have  not  yet  been  implemented  numerically. 
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FIG.  8.  Cross  sections  of  real  and  imaginary  parts  of  test  object  reconstruc¬ 
tions  using  data  from  five  wave  numbers.  Linear  reconstructions  were  ob¬ 
tained  using  retransmitted  fields  in  a  homogeneous  medium,  while  nonlinear 
reconstructions  were  obtained  using  retransmitted  fields  in  a  homogeneous 
medium  and  in  a  background  cylinder  of  compressibility  contrast  0.01.  (a) 
Cross  sections  of  real  parts,  (b)  Cross  sections  of  imaginary  parts. 


The  general,  nonlinear  version  of  the  eigenfunction 
imaging  approach,  as  defined  in  Eqs.  (24)  and  (44)  and  illus¬ 
trated  in  Fig.  8,  has  obvious  application  to  iterative  recon¬ 
struction  of  unknown  inhomogeneities.  In  such  reconstruc¬ 
tions,  the  total  pressure  field  can  be  estimated  at  each 
iteration  from  a  numerical  solution  of  the  direct  scattering 
problem  for  the  currently  estimated  inhomogeneity,  and  an 
eigenfunction  inversion  can  be  performed  using  this  pressure 
field.  This  procedure  can  then  be  repeated  to  obtain  im¬ 
proved  estimates  of  the  scattering  potential  until  convergence 
is  achieved. 

V.  CONCLUSION 

A  method  for  focusing  and  imaging  using  scattered 
fields  has  been  presented.  The  method  outlined  here  makes 
use  of  the  physical  properties  of  scattering  operators  by  using 
their  eigenfunctions  as  incident-wave  patterns.  The  eigen¬ 
functions  have  been  shown  to  provide  optimal  focusing  on 
pointlike  and  distributed  scattering  objects. 

The  inverse  scattering  scheme  presented  exploits  the  fo¬ 
cusing  properties  of  eigenfunctions  as  well  as  recent  analytic 
results  to  obtain  a  robust  and  efficient  means  of  quantita¬ 
tively  reconstructing  unknown  scattering  media.  This  new 
method  has  a  complexity  dependent  only  on  the  size  and 
complexity  of  the  scattering  medium.  Particular  cases  of  the 
method  provide  improved  implementations  of  eigenfunction 
focusing  and  filtered  backpropagation.  The  method  can  be 
implemented  for  any  background  medium  for  which  the  total 
acoustic  pressure  field  can  be  estimated. 

Another  particular  case  of  the  presented  method  results 
in  a  nonlinear  inverse  scattering  formula  that  yields  a  solu¬ 
tion  for  the  scattering  potential  q  in  terms  of  retransmitted 
fields  of  eigenfunctions  in  the  scattering  medium  and  in  the 
background  medium.  This  formula  has  been  demonstrated  to 
yield  improvement  in  accuracy  and  resolution  over  Bom  in¬ 
version. 
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The  ideas  reported  here  are  expected  to  be  useful  in 
further  studies  of  inverse  scattering,  adaptive  focusing,  and 
ultrasonic  imaging.  The  eigenfunctions  of  the  scattering  op¬ 
erator  A ,  whether  determined  by  iterative  retransmission  or 
by  numerical  diagonalization,  may  be  used  to  focus  through 
inhomogeneous  media  and  to  determine  optimal  incident- 
wave  patterns  for  inverse  scattering  experiments.  Also,  the 
products  of  fields  of  eigenfunctions  are  expected  to  form  a 
useful  basis  for  expansion  of  unknown  scattering  media  in 
iterative  reconstruction  algorithms. 
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